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Preface 


In 1979 a remarkable symposium "Glacier Beds: the Ice-Rock Interface" has been 
organized by the Canadian National Research Council. A wealth of excellent papers has 
been presented and discussed. Nevertheless, the final goal of formulating a general 
glacier sliding law seemingly had moved further away than before. Several speakers 
stressed this point; for instance Professor Hans Weertman, with his contribution "The 
unsolved general glacier sliding problem" or Professor Louis Lliboutry with the 
statement: "The best model is that which fits the largest number of observations and has 
the least complexity". Since this symposium, the interest in the sliding problem has 
grown steadily. It was roused especially by spectacular natural events: glacier surges, 
sudden slides of hanging glaciers, the disintegration and drastic retreat of tidewater 
glaciers and the "fast flow" of certain polar ice streams. Progress in understanding glacier 
sliding has been achieved on theoretical grounds and by applying highly developed 
numerical methods and experimental techniques. While no really "general" sliding law is 
yet in sight, several aspects of sliding are being better understood, for instance the 
influence of bed separation on sliding over bedrock. 


In the present study, "Glacier sliding over sinusoidal bed and the characteristics of 
creeping flow over bedrock undulations" the ice flow close to bedrock undulations is 
investigated in detail. The sinusoidal bed is the one with least complexity, and yet, it 
comprises a good model for typical glacier-polished bedrock, where small roughness 
elements are absent. Dr. Hilmar Gudmundsson derives the velocity field and stress field 
analytically for a linear-viscous medium and small bed roughness. He then extends the 
study to include large bed roughness and non-linear rheology by applying a suitable 
numerical method. He establishes the conditions under which local extrusion flow takes 
place. Beyond a certain limiting value of bed roughness a new phenomenon appears: 
circulating flow in the troughs of the sinusoidal bed. In this case the ice can no longer 
escape from the troughs. 


Besides giving a full account of the flow field near the interface of sliding, this study may 
thus help to understand certain unexpected disturbances such as local extrusion flow 
which have been encountered near the bed in measurements of borehole tilt. Furthermore, 
it will permit to assert the conditions under which ancient ice can be preserved in bedrock 
troughs. At a different scale, it shows the typical distribution of strain rates in glaciers 
which slide through overdeepenings. This last issue supplements the results of a 
subsequent study by Dr. Gudmundsson (No 131 of this publication series) of ice flow at 
a glacier confluence. 


This research has been funded by the Swiss National Science Foundation. 


Almut Iken 
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Abstract 


The characteristics of creeping flow close to bed undulations, and the form of the 
sliding law, in the absence of friction, bed separation and regelation, are analyzed. 
This is done by an analytical and numerical treatment of a non-linear medium 
flowing over a sinusoidal bed. The main theoretical results of relevance obtained 
to date, which are usually only valid in the case of a linear rheology, are given, 
and emphasis is placed on what can be learned from them about flow behavior 
close to the bed. An estimate of the pressure variation along the bed and of an 
effective viscosity is used to obtain a sliding law for a non-linear medium valid 
in the limit of small roughness and infinitely thick glaciers. Solutions based on a 
second order perturbation analysis for the velocity field are presented and examined 
for a sinusoidal bed. It is found that close to bedrock undulations, depending on 
the amplitude-to-wavelength ratio, two different regions of extrusion flow thay arise. 
Above the crest of the sine wave a region of local maximum, and within and above 
the trough a region of local minimum of the vertical velocity can develop, and exact 
criteria for the appearance and disappearance of these stationary points are given. 
Extrusion flow will cause a reversal of bore-hole inclination profiles close to the 
bedrock. This has been observed in nature but its cause has not so far been fully 
understood. 


Numerical calculations are performed to extend the analytical results to the case of 
non-linear rheology and a strongly undulating bed. The general form of the sliding 
law for a sinusoidal bed for every possible roughness and n value (n is a parameter 
in Glen’s flow law) is found, and the effect of the ratio of glacier thickness to bedrock 
wavelength is analysed. Extrusion flow is found to become increasingly important 
as the flow gets progressively more non-linear. For high roughness values a flow 
separation occurs, z.e. the main flow sets up a secondary flow circulation within the 
trough, and the ice participating in this circular motion theoretically never leaves 
it. 
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Zusammenfassung 


Die FlicBeigenschaften eines kriechenden Mediums in der Nahe von Bettuneben- 
heiten mit analytischen sowie mit numerischen Methoden untersucht. Dabei wird 
angenommen, da zwischen dem Bett und dem Medium keine Reibung stattfindet. 
Nur der Fall indem sich das Medium vom Bett nicht ablést, wird betrachtet. Regela- 
tion wird ignoriert. Ein Uberblick iiber die bisherigen theoretischen Ergebnisse, die 
sich in der Regel auf den linearen Fall beschranken, wird gegeben. Eine Abschatzung 
der Druckschwankung am Bett und der effektiven Viskositat wird benutzt um zu 
einem Gleitgesetz zu kommen das fiir den nicht-linearen Fall und kleine Rauhig- 
keiten giiltig ist. Neue Lésungen fiir das Flie8- und das Spannungsfeld, die auf 
einer Stérungstheorie zweiter Ordnung basieren, werden prdsentiert, und die wer- 
den fiir den Fall eines sinusférmigen Gletscherbettes untersucht. Es zeigt sich, daB 
je nach dem Verhdltnis der Amplitude der Bettunebenheiten zu deren Wellenlange, 
sich bis zu zwei Zonen bilden kénnen, an denen die horizontale Geschwindigkeit mit 
der Tiefe zunimmt. Oberhalb des héchsten Punktes des sinusférmigen Bettes kann 
ein lokales Geschwindigkeitsmaximum entstehen, und oberhalb des tiefsten Punktes 
kann sich ein lokales Geschwindigkeitsminimum bilden. Genaue Kriterien fiir dieses 
FlieBverhalten werden gegeben. Diese Art vom Flie8verhalten ist schon in der Natur 
beobachtet worden, aber bisher waren keine theoretischen Erklarungen bekannt. 


Anhand von numerischen Berechnungen werden die analytischen Resultate auf den 
nicht-linearen Fall und groBe Rauhigkeiten erweitert. Die allgemeine Form des Gleit- 
gesetzes fiir ein sinusférmiges Bett in Abhdngigkeit des Rauhigkeitsparameters und 
dem Grad der Nichtlinearitat wird bestimmt. Ebenfalls wird der Einflu8 einer end- 
lichen Gletschermachtigkeit auf das Gleiten untersucht und quantifiziert. Fiir sehr 
hohe Rauhigkeiten kann in einer Ubertiefung eine zirkulierende Strémung einset- 
zen. Das Kis, das sich in einer solchen zirkulierenden Strémung befindet, wird die 
Ubertiefung nie verlassen. 
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CHAPTER 1 
Introduction 


1.1 Basal sliding and flow close to bedrock 


Glacier flow consists of three elements: internal deformation of the ice, the rheology 
of may be described mathematically by Glen’s flow law, basal sliding, and deforma- 
tion of the underlying substrate. The sliding motion is often responsible for a large 
part of the flow velocity of temperate glaciers, and is therefore a subject of vital 
importance to the understanding of glacier behavior. The dependency of the slid- 
ing velocity on physical conditions at the ice-bed interface is still an open question. 
There have been many proposals for sliding laws, but so far there seems to be no 
general consensus on what are the pertinent variables governing the sliding velocity, 
let alone the form of the sliding law. Without a theoretically and experimentally 
well-founded sliding law, realistic basal boundary conditions for a model of glacier 
flow are not known and have to be introduced in an ad hoc manner. 


For a numerical flow model, sliding is a sub-grid process (for the grid densities which 
are currently practicable); sliding law is then expected to give the contribution of 
all sub-grid processes to an average basal velocity u, (where the averaging distance 
is of the same order as the distance between two grid points), as a function of some 
quantity calculated by the numerical model, as an average over a characteristic grid 
dimension such as the average basal shear stress 7. In general the sliding law has 
therefore the form: 


T. = f (us) (1.1) 


where f is some unknown functional representing the sliding law. The search for a 
sliding law is a search for a matching condition for the local (sub-grid scale) flow 
and the large-scale flow as has been stressed by Fowler (1986). Properties of the 
local flow determine and give the boundary conditions, i.e. the sliding law, for the 
large-scale flow. 


For any particular glacier the physical conditions at the bed-rock interface as well 
as factors known to influence the sliding velocity ——- the form of the bed, water 
pressure, debris concentration etc. —— are in general, not known. It will hence be 
difficult to obtain information on the functional dependency of the sliding velocity 
on for example the bed roughness spectrum through measurements. This illustrates 
the usefulness, indeed the need for a theoretical work on sliding. 


A short overview of previous work on the form of the sliding law and on flow char- 
acteristics close to bedrock undulation, in the presence of sliding, is given in Sec- 
tions 2.2 and 2.3 respectively. 
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The only process so far analysed in such detail that it can be considered to be 
formally understood is the sliding of a Newtonian fluid over a perfectly smooth 
wavy surface in the limit of small roughness, which can be solved with standard 
perturbation methods. The sliding velocity for high roughness values, even for this 
highly idealized model, is not known. It is not even known what roughness values 
should be considered to be high, i.e. how accurate the perturbation solutions are 
as a function of r, where r is the amplitude/wavelength ratio. No bounds on the 
accuracy of linear, as well as non-linear sliding laws using Glen’s flow law, proposed 
so far, are known and it is possible (and will be shown in Sec. 6.6.1 to be true) that 
the range of r values where the sliding laws can be used accurately decreases as the 
non-linearity of the viscous medium increases. 


Sliding over hard beds, as opposed to sliding over sediments or soft beds, is an 
important sliding mechanism in the Alps. Frictionless sliding over a sinusoidal bed 
without bed separation is possibly the conceptually simplest model of sliding that 
can be thought of. It is difficult to see how progress toward a “general” sliding 
law for hard beds, including the effects of friction and bed separation, can be made 
before a proper theoretical understanding of the more simple type of model has been 
obtained. A numerical approach seems to be the only possibility of addressing this 
problem in its full generality, since horizontal as well as vertical stress gradients 
are expected to become important close to bedrock protubances and the full set of 
equations has to be solved. Three dimensionless quantities characterise the problem 
(cf. Sec. 3.1); two geometrical numbers describing the roughness of the bed and 
the relative thickness of the glacier, and one number describing the degree of the 
rheological non-linearity. A three-dimensional parameter space must therefore be 
investigated, involving repeated solution of a non-linear equation system, so that 
considerable computer resources are needed. This is most probably the reason why 
this problem has not been tackled in its full generality by numerical means so far, 
although valuable insights have been gained from numerical work done for several 
isolated cases (Raymond, 1978; Meyssonnier, 1989; Schweizer, 1989). This numerical 
work will be discussed in Sec. 2.2 and Sec. 2.4. 


Most flow solutions, including solutions for large-scale datum flow (typical wave- 
length of bedrock undulation large compared to ice thickness) as well as short-scale 
solutions (typical wavelength of bedrock undulation comparable to ice thickness) 
do not calculate sliding from first principles but regard it as resulting from the 
local-scale flow (typical wavelength of bedrock undulation much smaller than ice 
thickness) in some unknown manner, and circumvent the problem by using an ad 
hoc sliding law usually of Weertman type 


up = Cr,” (1.2) 


where m and C are adjustable parameters. Although this law has mathematically 
the same form as Glen’s flow law it does not have the same experimental status 
(Lliboutry, 1968). The important point here is that the amount. of basal sliding is 
not an output but an input of most flow solutions. This makes perfect sense for 
large-scale flow models but may be somewhat. questionable for short-scale solutions. 
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Given the current state of affairs there seems, however, to be no other way of tackling 
the effect of sliding on the flow. 


Little is known about the local flow properties close to bedrock undulations in the 
presence of sliding even for the linear case, and almost nothing for a non-linear 
medium (cf. Sec. 2.3). Because of the inherent complexity of the problem, only a 
few analytical solutions exist. These solutions often apply to somewhat idealized 
conditions at the rock bed but nevertheless give a valuable insight in to the nature 
of the flow. Again only flow properties for small roughness are known and numerical 
work has been limited to a few cases. 


Measurements of the deformation of bore hole give information on the rheological 
properties of ice. Since shearing is usually concentrated at the lower most section, 
directly above the bed, it is important for the interpretation of bore-hole deforma- 
tion data to know how ice deforms around bedrock protrubances. Some surprising 
observations, which require an explanation have been made. Extrusion flow close to 
the bed has been demonstrated by bore-hole measurements (2.e. close to bed the ve- 
locity increased with depth) (Hooke et al., 1987), which could possibly be caused by 
the presence of bedrock undulations. Extrusive flow has also been observed within 
subglacial sediments (Blake et al., 1992). 


Flow through an overdeepening has been the object of increased interest with the 
recent finding in the Alps of a corpse buried in ice, which was dated to be about 
5000 years old. The fact that the corpse was found to be lying in an overdeepening 
suggests a way of how it could have been situated there for such a long time. In 
this context it is especially interesting to know to what extent ice flows over an 
overdeepening without activating the ice within it. This will of course be dependent 
on the depth-to-width ratio. A sinusoidal bed is a convenient idealization of this 
type of geometry. To what extent the ice lying in the overdeepening participates in 
the general flow and for how long it stays within it depends on the exact boundary 
conditions, but a lower limit to the time it stays there and an upper limit on the 
flow velocity is obtained by assuming no friction at all, t.e. perfect sliding. 


1.2 Goals of the study 


The goals of this study are twofold: (i) to investigate and determine the dependency 
of the sliding velocity on bedrock roughness and rheological non-linearity, and (ii) 
to analyse the characteristics of glacier flow close to bedrock undulations in the 
presence of sliding. 


To this end a study of an idealized glacier having a sinusoida] bed given by 29 = 
asin kx, where zo is the vertical position of the bed, a the amplitude and & the 
wavenumber, will be done. The single wave roughness is defined as the ratio of 
the amplitude to the wavelength and will be denoted by r, ze. r := a/A, where 
A = 27/k. Another convenient measure of the roughness is the local bed-slope 
parameter € defined as ¢ := ak. The ratio of the glacier thickness h to the bed 
wavelength is expected to influence the sliding behavior. The thinness parameter 6 
is defined as 6 := (kh)7!. 
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Attempts should be made to find answers to questions such as: 


e Is wy x 1/e"*! for ¢ > 0, and if that is true, for what range of € values can 
this relation be used? Does the range depend critically on n? 


e What are the effects of 6 on us other that those entering the sliding law through 
7? Does 6 influence uy significantly or can it for practical purposes be ignored? 


e How does uy, depend on « for ¢ ¢ 1? 
e Is extrusion flow close to bedrock undulations possible? 


e How does ice flow through an overdeepening? At what € values does ice 
effectively remain within the trough? 


e Is flow separation possible? 


Another goal of the study is to test the feasibility of the FE method for flow calcu- 
lations involving sliding, and to develop programs that automate the calculation of 
sliding velocities for a general bed geometry. 


1.3. Organization 


In Chapter 2 a specification of the problem that forms the subject of this work is 


given. Previous work is discussed, with an emphasis on the linear first order theory 
of Nye and Kamb. 


Chapter 3 focuses on what is known about the form of the sliding law for a non- 
linear medium, and what can be learned from simple dimensional arguments. An 
estimate of the effective viscosity is used to show how strongly the sliding velocity 
depends on the slope parameter e. . 


The second order solution of Morland for the flow along the sole-bed interface. is 
extended in Chapter 4 to the region above the bed-line. It is shown that the solution 
so found displays a number of interesting properties. 


For a non-linear flow law of the type mostly used in flow modeling, i.e. Glen’s flow 
law, analytical solutions have not been found. The only way at present of obtaining 
a better understanding of this important problem seems to be through the use of 
numerical models. In Chapter 5 a testing of the correctness of numerical results 
obtained with the FE program MARC is described. This program was used for 
numerical calculations of flow over a perfectly lubricated sinusoidal bed. 


In Chapter 6 the numerical results are presented and discussed. 


A summary of the results is given in Chapter 7. 
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CHAPTER 2 


Previous Work 


In this chapter an exact description of the idealized problem of a glacier flowing 
over a perfectly lubricated sinusoidal bed is given. This includes a mathematical 
description of ice as a highly viscous isotropic and homogeneous material. Previous 
work on sliding, with emphasis on sliding without bed separation and regelation, is 
summarized. 


2.1 Description and specification of the problem 


2.1.1 Ice dynamics 


Laboratory experiments and field observations suggest a constitutive law for ice of 
the form 


. + (n-1)/2 + 
6 = Aor, ij (2.1) 
. . ‘ . . 
where €;; are strain rates, 95; deviatoric stresses, 
’ 1 
9%; = O73 — 3 fii Tke (2.2) 
: . . . . 
o,, the second invariant of the deviatoric stress tensor 
# 1 , t 
O11 = 595%; (2.3) 


and 4;; is the Kronecker delta. A and n are parameters determined by measure- 
ments. Unfortunately the values of A and n are only approximately known. Often 
n is assumed to have the value 3, although the experimental basis for this is not 
strong. A, which is sometimes called the softness parameter, is a thermodynamic 
property and varies strongly with temperature. Other factors such as ice fabric, ice 
composition and chemical state and, for temperate ice, moisture content, are also 
known to be important and to influence the value of .4, e.g. Lile (1978), Budd and 
Jacka (1989), van der Veen and Whillans (1990). Sometimes the hardness parameter 
B = A~'/" is used instead of A. This type of constitutive law is known in metallurgy 
as Norton-Hoff’s power law, in general fluid dynamics literature as the Ostwald-de 
Waerde or Reiner power law, but as Glen's flow law in the glaciological literature 
(Hutter, 1991). Glen (1952, 1955) and Steinemann (1958a, 1958b) showed for the 
first time the necessity of a non-linear flow law for ice. Generalisation of Glen's 
experimental results — that were done with a stress state of uniaxial compression 
— by Nye (1953, 1957) resulted in the above shown flow law. 
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All glacier flow must obey the fundamental laws of conservation for a physical sys- 
tem, 2.e. the conservation of mass, linear and angular momentum, and of energy. 
Since the density of ice can be assumed to be constant, conservation of mass is 
equivalent to conservation of volume. 


Continuity equation (mass conservation): 


Ui = 0. (2.4) 


Equation of motion (conservation of linear momentum): 


0:35 + PGi = 0. (2.5) 
The acceleration term is ignored since it is truly small compared to other terms. 


Energy equation (conservation of energy): 


DT 
PORE = KT x + Oi; ij, (2.6) 


where D/Dt is the substantial derivative!, c the heat capacity, and K the thermal 
conductivity. In Eq. (2.6) use has been made of one constitutive equation, i.e. 
Fourier’s law of heat conduction. Since the ice will be assumed to be temperate, the 
temperature will be determined by the Clausius-Clapeyron relation, and Eq. (2.6) is 
inappropriate. The energy released by viscous dissipation will increase the moisture 
content and not raise the temperature. The softness parameter A is expected to 
be a function of moisture content (Lliboutry, 1976) but this will be ignored here. 
Only spatial variations of A within the basal layer will affect the results presented. 
It is not clear that the influence of spatial variation of moisture content on A will 
dominate other factors on which A depends, such as salt, air-bubble and debris 
content, which can either increase or decrease the value of A, and which will also be 
neglected. A theoretical formulation of temperated ice masses that accounts for the 
interaction of ice and water has been given by Hutter (1993) and numerical work 
has been done by Hutter et al. (1988). 


2.1.2 Problem definition 


Consider a glacier, thought of as a highly viscous material, sliding slowly and steadily 
over a sinusoidal bed. Inertia forces are negligible. The sliding is perfect in the sense 
that there is no tangential traction at the bed. The beds exert only normal forces. 
No bed separation occurs and regelation is neglected.?, The problem is depicted in 
Fig. 2.1. The bed profile is given by 


z= 2 = asinkz. (2.7) 





'D/Dt := 0/dt+v;0/r;. Other names for the substantial derivative are the particle derivative 
and the material derivative. 


?For real glaciers regelation is expected to be important at wavelengths smaller than the tran- 
sition wavelength (Weertman, 1957; Weertman, 1964; Weertman, 1979). 
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Fig. 2.1: Flow over a sinusoidal bed. The coordinate system makes the angle a with respect 
to the horizontal. The vertical position of the bed line zp, ts given by: z, = asinkx. The sine 
wave has the wavelength A = 27/k, and amplitude a. The surface velocity is denoted by u, 
and the sliding velocity by uy. / is the glacier thickness. 


2.2 Previous work on the form of the sliding law 


Weertman (1957) was the first to put forward a theory of glacier sliding. Later 
he refined and defended his model (Weertman, 1964; Weertman, 1971). He distin- 
guished between two fundamentally different physical processes. The first is rege- 
lation, whereby pressure variations around a bed obstacle cause a melting on the 
up-stream side and refreezing at. the down-stream side. The second is local viscous 
deformation of ice around bed obstacles.’ For a given driving stress 7 Weertman 
found the regelative and the viscous components of the ice velocity to be propor- 
tional to ™% and 7,’ respectively, where n is the exponent in Glen’s flow law. 


Lliboutry (1959, 1968, 1979, 1987b) argued that ice deformation, where the ice sepa- 
rates from the bed, forming water filled cavities, is an important sliding mechanisin. 
Bed separation reduces friction and causes larger sliding velocitics.’. Variations of 
water pressure cause a change in the extent of bed separation and must be linked 
to changes in sliding velocity (Uken, 1981). Measurements have confirmed this (ken 
and Bindschadler, 1986). The magnitude of sliding with bed separation must de- 
pend on the water pressure in one way or another. Theoretical work has been done 
by Fowler (1986, 1987) that. produced a complex sliding law for bed with power- 
law roughness spectra r(A) o AY. Further work has been done by Kamb (1987), 





*SWeertiman used the term enhaneed plastic flow to describe this mechanism that T call local 
mscous deformation. I do not use the word plastic because the flow is uot plastic (for 2 > 1 Glen's 
(low law represents a pscudo-plastic stress-strain relationship (White, 1991, p. 24-25)), and the 
word enhanced does not seem necessary since it is not clear with respect to what. the flow is being 
enhanced. 

‘Touse the term bed separation instead of the word cavitation since cavitation is used in Muid 
mechanies to denote another physical mechanics, namely the spontancous creation of air bubbles 
within a fluid as the {nid pressure becomes locally lower than the vapor pressure (Hutter, 1983). 
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Humphrey (1987) and Schweizer and Iken (1992). A recent overview has been given 
by Kamb (1993). 


Various attempts to estimate the relative importance of regelation and ice defor- 
mation with and without bed separation have all shown that regelation is only 
expected to be of importance at small distance scales (Weertman, 1957; Nye, 1969; 
Fowler, 1986). For a given ratio of the amplitude of the bedrock undulations to the 
distance between them, sliding due to regelation is inversely proportional to wave- 
length, but sliding caused by ice deformation increases linearly with wavelength. 
This is the motivation for Weertman’s notion of a controlling obstacle length, where 
the contribution to sliding by regelation and ice deformation is the same. By sim- 
ple dimensional analysis and assuming that the contribution of regelation and ice 
deformation (without bed separation) can be added to the sliding velocity, one finds 


Up p A Cok 


ts = a + “ly? (2.8) 
where 4 is the wavelength of the undulation (or the distance between the obstacles), 
7 is the viscosity of the ice, which for a non-linear material could be substituted by 
an effective viscosity, Cp is the Clausius-Clapeyron constant, K the mean thermal 
conductivity of ice and bedrock, and L is the latent heat of fusion per unit volume 
of ice. Gj and ¢; are some unknown constants, that in general will be some unknown 
function of all the dimensionless numbers and ratios entering the problem, such as 
the roughness r := a/A, the ratio of the wavelength to the thickness of the glacier, 
and a possible non-dimensional number entering the flow law. 


Nye (1969, 1970) and Kamb (1970) considered the flow of a linear viscous medium 
over a perfectly smooth sinusoidal bed without bed separation, and derived a first 
order perturbation solution for the flow and for the stress field. This solution will 
be discussed in Section 2.3.1. For a single sine wave Nye and Kamb found 


Up _ 1 2Co(K, + Kg)k 
T nek Le? ; 
where k is the wave number (k = 27/X), € := ak is the local bed-slope parameter (€ = 
2mr), and K, and Kg are the thermal conductivities of ice and bedrock respectively. 
The first term on the right-hand side represents sliding due to pure viscous flow and 


the second sliding caused by regelation. The equation can only be used for e¢ < 1. 
The transition wavelength A, is given by 


(2.9) 


ie 822Co( Ky + Kg) 
i L 

and is of the order of 0.5m. Morland’s (1976a) second order linear perturbation 
theory gave the same result. If roughness is absent at wavelengths smaller than 
X,, Tegelation can be ignored in the calculation of sliding velocities, and the sliding 
velocity is proportional to 7;. Conversely, if the bed roughness is concentrated 
at short wavelengths (A < 4,), then up «x 7. For a bed consisting of a single 
sine wave with \ = \,, Weertman (1957) showed that wy o Tee: The heat 
transfer associated with the regelation process was calculated by assuming it to be 


(2.10) 
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caused by pure conduction through rock and ice only, which can be justified to some 
extent by the thinness of the water film estimated to be about 1 zm (Nye, 1967; 
Lliboutry, 1968; Fowler, 1981). A more accurate estimate of melting and refreezing 
by regelation is difficult as the ice salt and air-bubble content, water inclusions, 
moisture content and surface properties of the bedrock must be accounted for (Drake 
and Shreve, 1973; Chadbourne et al., 1975; Morris, 1976; Shreve, 1984; Meyssonnier, 
1989). 


Kamb (1970) developed an approximate non-linear theory of sliding, which will be 
discussed briefly in Sec. 3.3. He made the assumption that the effective viscosity is 
a function of the distance above the bed only, and found that the transition length 
is somewhat smaller for n > 1 than for n = 1, where n is a constant in Glen's flow 
law. 


Fowler (1981) made a rational dimensional analysis of the sliding of glaciers in 
the absence of bed separation, including the non-linearity of the flow law. For 
the transition length he found a value of the order of 0.01m and concluded that 
the regelative component of the ice velocity may therefore be neglected, provided 
roughness is absent at lowest wavelengths. 


By using a variational theorem for non-Newtonian flow found by Johnson (1960), 
Fowler (1981) derived a sliding law valid for a non-linear medium, concluding that 
the sliding velocity u, should be proportional to the ratio 6/e"*', where 6 = 1/kh 
and h is the glacier thickness. A sliding law with the same functional dependence 
on € had earlier been derived by Lliboutry (1968) using simple physical arguments. 


That uy is so highly dependent on « for n = 3, the value of n mostly used for ice 
flow modeling, is striking. This has, however, been questioned by Schweizer (1989) 
who, through FE calculations, found up to be less dependent on € than up « 1/e"*! 
indicates. 


Raymond (1978), in what seems to be the first numerical calculation of this prob- 
lem, calculated the sliding velocity for a perfectly lubricated sinusoidal bed with no 
regelation using Glen’s flow law. His calculations were limited to a few ¢ values and 
they do not give the sliding velocity as a general function of ¢€, 6, and n, but they 
stand in an agreement with up, « 1/e"*?. 


Meyssonnier (1983) did a comprehensive numerical study of this problem, and com- 
pared numerical calculated sliding velocities with analytical estimates. His impor- 
tant work will be discussed in Section 3.2. 


2.3. Previous work on flow characteristics close to 
bedrock undulations 
2.3.1 The linear first order theory of Nye and Kamb 


Nye and Kamb, and later, in a somewhat more general way, Morland (1976), have 
found an approximate solution to the problem of a highly viscous Newtonian fluid 
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sliding over a smooth bed. The assumption of a Newtonian fluid was necessary in 
order to be able to use standard perturbation methods. Morland has also considered 
the effect of friction on the flow (Morland, 1976b). 


The method of Nye and Kamb will now be briefly outlined, mainly in order to point 
out the approximations that had to be made in order to obtain the solution. 


A solution for the velocity v = (u, w) and the pressure p is sought in the form 


u= up + > e'u,(z, z), (2.11) 
1=1 
w= do e'ui(z, z), (2.12) 
i=1 
p= Ye'p(a,2), (2.13) 
t=1 


where € = ak. One may think of € as some small number, ¢ < 1. The bed profile is 
given by 
z= 2% := asinkr =: ef(z). (2.14) 


At the upper boundary where z = z, with z,; > 4, a shear stress 7 is applied. 
This is the driving force of the motion, not gravity, which is not present. The idea 
here is that the force of gravity acting in the layer close to the bed will be negligible 
compared to the force that the overlaying ice is exerting on it; a good approximation 
for a thick glacier, that is if 6 ~ 0. This means that close to the upper boundary 
the horizontal velocity varies linearly with z. On the lower boundary, zp = asin kz, 
there are two conditions imposed: the velocity normal to the bed is zero,” and there 
is no tangential traction. This means that 


d 
Sig OA yelp , on 2 2h, (2.15) 
dz 
and that ; 
n= 5 (Oe —0,,)tan2B(r) , on = 29: (2.16) 
where tan G(r) := dzo(xz)/dx. These exact boundary conditions are now applied 


not at the real boundary z = zo, but at z = 0. By inserting zy) = ef(x) and the 


expressions (2.11) and (2.12) into Eq. (2.15), and by only taking terms up to order 
€, one obtains, since 2 is O(e): 


—usf (x) + w(x, 0) = 0. (2.17) 


The second boundary condition can be approximated in a similar way. Since both 


’ 


Onz = (Oxz ~ 0:z)/2 and tan 2G(x) are O(c), the boundary condition (2.16) becomes 
a(x, 0) = 0; (2.18) 


writing o;, as 
Oz2 = €0,; + Ole”). (2.19) 





*In the presence of regelation this is not correct and there is a non-vanishing velocity component 
normal to the bed. Regelation was taken into account by Nye and Kamb, but it is ignored here. 
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Or if expressed in terms of v = (u, w) 


f oO “ 
w, = uf (z) and =-wuf (zr), on z=0. (2.20) 


The expressions for the velocity and the pressure are now substituted into 
nV’v =Vp and V:v=0 (2.21) 
the conservation laws of linear momentum and mass. 


After having collected the first order terms the solution is obtained by usual Fourier 
integral methods. In the special case of a bed consisting of a single sine wave these 
first order perturbation solutions are: 


To 


Uy = ae (2.22) 

Uz(Z, Zz) = Uy + uak*ze~™ sin kz + O(e?), (2.23) 
v,(z,z) = upk(1 + kz)e~*acos kz + O(e?), (2.24) 
P(Z, Z) = Poo + 2nusk?e~**a cos kz + O(e?), (2.25) 
0,,(Z,2) = —0,, = 2nupk?ze-**acos kx + O(e?), (2.26) 
Ox2(Z, Z) = —2nupk?ze~*asin kx + O(e?), (2.27) 
T(x, z) = 2nuak*ze-** = QT» —~e** + O(e?), (2.28) 


a 
where 7 := Vom which is sometimes called the effective stress, and p. is the 
pressure applied at the upper boundary of the medium. 


These expressions have some interesting features. One of them is the fact that the 
second invariant of the deviatoric stress tensor, Eq. (2.28), shows no dependence on 
z. This will of course also apply to the second strain-rate tensor invariant. 


Another interesting feature of the linear solutions above is the occurrence of extru- 
sion flow, there is a region close to the bed where the horizontal velocity increases 
with depth. At the point ka = 37/2+ 27! where ! € N and z = zmin := 1/k, vz has 
a local minimum: 


Umin = Us(Z = 3n/(2k), z = 1/k) = ug(1 — ak/e). (2.29) 


From here downwards to the bed the horizontal velocity increases. This is a re- 
markable fact, which deserves a closer inspection. Notice that since ak < 1, it 
follows that zmin >> a. Extrusion flow, a term introduced by Demorest (1941, 1942) 
has been a subject of some debate in the glaciological literature. On theoretical 
grounds it can easily be shown that a global extrusion flow, that is an increase of the 
horizontal velocity with depth throughout an entire glacier, is impossible since the 
overlying mass will then experience a force in the main direction of flow, which is not 
counterbalanced by any other force, leading to an accelerating velocity (Nye, 1952). 
There are, on the other hand, claims of extrusion flow having been directly observed 
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by bore-hole deformation measurements (Hooke e¢ al., 1987), and by observations 
within sub-glacial caves close to the bed-rock interface (Carol, 1947). Arguments 
supporting (global) extrusion flow based on mass-balance measurements have also 
been given (Streiff-Becker, 1938; Seligman, 1947). These arguments must, however, 
be considered to be rather weak, because no direct information on the velocity field 
was available. 


Because the solution of Nye and Kamb is only approximate, one might be tempted 
to believe that this occurrence of extrusion flow may only be an artifact due to its 
inexact nature, and that a more complete treatment where second order perturbation 
terms are included will eliminate this feature. The inclusion of gravity as the driving 
force of the motion, instead of a given shear stress at the upper boundary, will modify 
the velocity profile so that it no longer varies linearly with depth at some distance 
above the bed and that would also affect this situation.® A somewhat more realistic 
theory is needed to clarify these matters. In Sec. 4.1 the second order solution 
of Morland (1976a) will be used to calculate the velocity and the stress field as 
functions of z and z (Morland only gave the solution along the bed, i.e. z = 0) and 
the solutions will be analysed with respect to the possibility of extrusion flow. 


2.4 Numerical work 


Raymond (1978), Meyssonnier (1983), and Schweizer (1989) did FE calculations of 
flow over a simple sinusoidal bed. 


Raymond found the excess of velocity in comparison with a no slip situation to 
depend on height above bed, which he ascribed to an effective softening of the basal 
ice for n > 1. He also found the deformation to be concentrated more closely near 
the bed for non-linear than linear behavior. 


Meyssonnier obtained a point of maximum relative horizontal velocity that was 
situated above the peak of the sine wave. This point did not appear in all his 
calculations, for which the reason is not clear. 





SIt is interesting that neither Nye nor Kamb mentioned the feature of there solutions which 
suggests the existence of extrusion flow for an infinitely wide glacier at small asperity and roughness 
values. My guess is that they thought of it as an artifact of no relevance for glacier How. 
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CHAPTER 3 


Sliding Law for a Non-Linear Medium 


No analytical solutions exist for the flow over a sinusoidal bed of a power law medium 
with n # 1. Despite this fact some information on the form of the sliding law can 
be obtained through the use of dimensional analysis and variational theorems. 


3.1 Dimensional arguments 


The relevant variables are: the glacier thickness (h), the wavelength of the sinusoidal 
bed (A), the amplitude of the sine function (a), the mean slope (a), the driving stress 
(7 = prghsina) — which is equal to the mean drag — and the viscosity 7. The 
sliding velocity (us) can only depend on these variables. That is 


Up = us(a, A,h, 7,7). (3.1) 


There are two dimensionless ratios that can be formed with these variables: the 
roughness, r := a/A and the asperity, ¢ := A/h, were h is the glacier thickness.' It is 
sometimes more convenient to use, instead of r and ¢, the parameters: ¢ := ak, called 
the local bed-slope parameter’, and 6 := 1/kh, which I call the thinness parameter. 
For an infinitely thick “glacier” the thinness parameter is equal to zero.’ 


The sliding velocity will in general depend on both parameters. If, however, the 
glacier thickness, h, is much larger than the wavelength \ and the amplitude a, 
it is physically reasonable to expect that the ice close to the ice-bed interface only 
responds to the applied stress caused by the overlying mass of ice and that it does not 


'Fowler (1979) introduced the word asperity for ¢ and used the word corrugation for ¢. To call 
the ratio of amplitude to wavelength roughness (sometimes multiplied by some numerical factor), 
and use the symbol r for it, has been done by many authors, e.g. (Nye, 1970; Kamb, 1970; 
Lliboutry, 1968; Schweizer, 1989; Schweizer and Iken, 1992). 

2¢ is the so called small slope parameter which Nye and Kamb introduced (Nye, 1970; Kamb, 
1970). Here I simply call it the slope parameter since in what follows it will not necessarily be 
small, or the local bed-slope parameter to make it clear that. it is the local slope of the bed which 
is being referred to. 

3Theories of large-scale datum flow, e.g. shallow ice approximation (cf. (Hutter, 1983)) often 
use the symbol 6 for the ratio of the vertical to the horizontal dimensions of an ice cap, which is 
then inversely proportional to the thinness parameter defined here and also denoted by 6. The 
motivation for defining 6 in this way here, is that the symbol 6 is usually thought to stand for 
some small dimensionless number, and in what follows 6 := 1/kh has to be small. The assumption 
of § < 1 is not a self-imposed restriction; it is difficult to see how sliding velocity, which must be 
some average property of the inner flow (i.e. the flow close to bed) can be defined if this does not 
hold. The assumption 6 < 1 can be seen as part of the definition of the problem. 
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sense the existence of a free surface, so that the thickness h will not appear explicitly 
in the sliding law. In this limit of an infinitely thick glacier w, is a nontrivial function 
of r only (nontrivial in the sense that this function cannot be determined by mere 
dimensional arguments) or equivalently of the slope parameter ¢€ := ak = 2Qrr, 
Dimensional arguments now give 


us = g(€)mA/N (3.2) 
or 


up = G(€)ma/n, (3.3) 


where g and g are some unknown functions and 7 is the average basal shear stress. 
Assuming that the rheological behavior of glacier ice can be approximated in suff- 
cient detail by Glen’s flow law 


» (n-1)/2 1 
ig = Ao O55 (3.4) 


/ . . . . 
where o,, denotes the second deviatoric stress invariant: 


’ 1 i i 1 f 
the sliding law becomes 
up = U,(e, 6, n) ATP A. (3.6) 


Notice that the rf dependence of uy follows entirely from dimensional arguments 
and that this is the only form possible if the effects of other physical processes such 
as regelation, are neglected. 


The unknown function U, can, and most probably will, depend on n and 4 since n 
and 6 are dimensionless. One can think of U, as a nondimensional sliding velocity. 


3.2. The sliding law in the limit as « — 0 and for 6 + 0 


It is possible to show how the sliding velocity, u,, depends on the (local bed) slope 
parameter ¢, and the parameter n in Glen’s flow law for 6 ~ 0 in the limit ¢ —> 0, 
if one assumes that the strain rates close to the bed are mainly determined by the 
kinematics of the problem, and that changing n will have much more effect on the 
stresses than on the strains. This will be done below. Somewhat similar arguments 
have been given by Kamb (1970) and Lliboutry (1968). The argument given suggests 
a certain form of the sliding law but does not prove that it must have this form. 
The result has already been proved by Fowler (1986, 1987) and does not, of course, 
have to be proved again. The arguments used by Fowler are, however, somewhat 
complicated, and it is instructive to see how the form of the sliding law, can be 
derived in a much simpler, albeit less rigorous, way. 


The idea is to express the pressure variation along the bed as a function of the 
amplitude, the wavelength, the sliding velocity and the viscosity. Since the rheology 
is non-linear the effective viscosity will be considered. It is assumed that ¢ <1. 
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The pressure along the bed‘ is then given by 
P(x) = pa + prghcosa + dp, (3.7) 


where 
Op(x) = cousnegke cos kx + O((ak)?). (3.8) 


Cy is some unknown constant, p, is the atmospheric pressure and mg is an effective 
viscosity that will be defined below. In order to see the plausibility of this expression 
notice that dp must vanish as € — 0. That dp(z) depends linearly on € can be thought 
to be the result of a Taylor expansion with respect to € where only the first term is 
retained. This will always give a correct result if ¢ is small enough.’ The product 
nuk gives the right dimension for the pressure. It would also have been possible 
to use the basal shear stress 7 to get the right dimensions, but then 7 would not 
have been a part of the expression. Notice also that k cannot be substituted by 1/a 
because if a + 0, dp must vanish (another way of seeing this is that if at» —a, dp 
must change sign). It is obvious that the variation of dp with z must be of the form 
cos kz if € is sufficiently small. 


Now the effective viscosity jeg will be estimated. 


As a column of ice, which extends from the bed up to some distance ! where the 
disturbance due to the bed undulations has disappeared, moves the distance Az, 
in the time At = Az/up, it is stretched/compressed to the length | + Az, where 
Az = Ardzg/dr. The (average) strain rate is given by 


1 Az 


<é,,>= Ard? (3.9) 


where the brackets are used to indicate that these are averaged strain rates over a 
given distance, that is 


: 1 /'. 
<6 = 7 i €4j dz. (3.10) 
0 
For the average vertical strain rate one gets 


1 Axcakcos kx 
Az/us 1/k 





< 6,2 >K + O((ak)”) = upak? cos kr + O((ak)’), (3.11) 


where use was made of the fact that | must be proportional to 1/k since the only 
parameters entering the problem that have the dimension length are a and 1/k, and 
that using a would again be incorrect since the strain rates must change sign as 
at+ —a.®° This expression could also have been obtained by calculating the average 
of the vertical strain rates over the whole glacier using the vertical strain rates from 


‘Note that the pressure at the bed is not equal to the negative of the local normal stress o,, 
i.e. p # —o,. This can be seen by writing the pressure as p = —$(on +o,), where gp is the local 
stress parallel to the bed, and noticing that o, — op = on - op # 0. 

5Qne could also do the following analysis to second order in e. I have done this but it did not 
give any new results. Note that the first order term is needed since Nye and Kamb’s solution shows 
that Op(z)/de # 0 for « = 0. 

®Notice that it is the slope of the sine curve, which has to be small if this analysis is to be 
correct. It is therefore « which has to have a value of less than unity, and not a/d = ¢/2r. 
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the theory of Nye and Kamb discussed in Sub-sec. 2.3.1. Doing this for the shear 
strain rates gives for <€z, > 


<€,, >= wyak’ sin kz. (3.12) 
The effective strain rate € is 
. 1. , : : 
= a fitis = e+e, (3.13) 


which now turns out to be 


é = uak* cos? kz + sin? kz = uyak?. (3.14) 


The Glen flow law can be written as 


; e(l-n)/n ; : 
07; = An “4 = 2Net Ei; (3.15) 
where 
e(l-n)/n 
Nef “= Sai” (3.16) 
Inserting (3.14) and (3.16) into (3.8) gives 
p(x) x |uyak?/A|'/" cos kz. (3.17) 


Force equilibrium requires 


sin @ 1 ‘ 
pigh aa) = ay 74 (3.18) 


where ¥ is a path along the sine curve, C'(y) is the path length, o is the stress tensor, 
and 7 is a unit normal vector 


Boks 1 —dz,/dx 
1+ (dz,/dx)? ( 1 , ey) 


A change of variables dr = ds(1+(dz,/dz)*)~'/? gives for the drag (the z-component 
of Eq. (3.18)) 


1 p dz 
Th = a (922 — Ose 7) de 
1 x dz ' 
LS oes pa hee . 3.20 
xd (o— + On2 — 0,, tan 3) dx (3.20) 


where tan 3 = dz,/dz. The last term on the right-hand side is of second order in €, 
and the exact boundary condition on zo (Eq. 2.16) shows that a,, is then of second 
order also. Up to first order in € the drag can therefore be calculated according to 


7 °= pighsina = ode Ce a (3.21) 
b= Pg me i P(Z, Zo Bs a : 
giving 
ck 2n/k , i 
Th = oo. (usak?/A)!/"ak cos? kx dz, (3.22) 
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where ¢ is an unknown constant. After integrating and rearranging terms one gets 


2Ar" j 
|(ak)"+*k| 


Up XX as (3.23) 


Expression (3.23) is a remarkable result. It brings out the asymptotic 1/e"*! be- 
havior of the sliding velocity as e — 0 and shows how strongly it depends on the 
bedrock roughness. It also encourages the definition of a new function s, which I 
call the sliding function 


et kus(e, 6, n) 


s(€,6,n) := DAT 


(3.24) 
One can think of the function s as the sliding velocity brought to a non-dimensional 
form, where the asymptotic behavior for ¢ — 0 has been accounted for. s is an even 
function of ¢. For every particular n, 6 + 0, and e < 1, s is a constant. In general, 
however, it is of course dependent on €. From (2.22) it is seen that s(0,0,1) = 1. 


Kamb (1970) found the same asymptotic behavior as in Eq. (3.23), but his derivation 
involved various assumptions which made it difficult to judge the general correctness 
of his findings. Kamb’s work will be discussed in Sec. 3.3. 


It wasn’t until Fowler (1979), by using a variational theorem (Johnson, 1960), came 
to this same conclusion, that Eq. (3.23) was given a sound mathematical basis. 


Lliboutry (1987b, p. 356) gave an argument somewhat similar to the one given 
above. His line of thought is however not easy to follow and one of his equations, 
which gives the basal shear stress (Eq. (13.39)) as a function of o,,, is incorrect. 


Schweizer (1989), in a numerical study, came to the conclusion that the sliding 
velocity does not vary as e~("+!), but somewhat more slowly. His results are also 
not in accordance with results based on a simple dimensional analysis (Eq. (3.6)). 
They must therefore be, to some extent, inaccurate. 


The function s(e, 6,n) can be expressed as a Taylor series with respect to the variable 
€ 
s(e,6,n) = y Cai(5,n) €**, (3.25) 
t=0 
where c,; are the Taylor coefficients. 


Fowler (1981), Lliboutry (1987a) and Meyssonnier (1983) were able to estimate 
co(0, 3). Meyssonnier found 


0.305 < co(0,3) < 0.338. (3.26) 
By using numerical methods to calculate the sliding velocity tor values of r in the 
range 0.01 to 0.05 and 6 in the range 0.25/m to 1/7, Meyssonnier was able to give 


an estimate of c2(0, 3) 
c9(0, 3) = 2.4. (3.27) 


He gave no estimate of the errors involved and did not perform a systematic inves- 
tigation of the 6 dependency of u,, but concluded that it is small. 
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Lliboutry (1993) used a variational method to obtain an upper bound on s for n = 3, 


His result is 
5(0,3) < 0.33839 + 3.688 e” + 0.169e° (3.28) 


Eq. (3.23) shows that s(e,6,n) as a function of n will be equal to cé™", where ¢ is 
some unknown number. Since s(0,0,1) = 1, é= é, giving 


s(,,n) = ee. (3.29) 


Nye’s and Kamb’s solution for the pressure p(z,z) shows that é = 2 for a linear 
Newtonian medium. Assuming that the effect of a non-linear flow law on the pressure 
variation along the bed can be described by a corresponding change in the effective 
viscosity, é in (3.8) will remain the same for different values of n. The estimate of 
the effective viscosity given above may, however, not accurately describe the effect 
that a change of n has on mq. The proper integration length / in (3.10) will for 
example depend on n if the deformation is concentrated in a different way near the 
bed for non-linear than for linear behavior. Important aspects of the non-linearity 
of the flow law such as these cannot be described simply by introducing an effective 
viscosity. 


3.3. Kamb’s non-linear sliding law 


The effective stress is, as was discussed on page 25, to first order in € independent 
of x. Kamb (1970) used this fact as a starting point for the development of a 
theory of sliding incorporating rheological non-linearity; assuming that since €;; is 
independent of z in the linear theory (where only first order perturbation terms are 
considered) it will also be approximately so for non-linear flow, in which case the 
effective viscosity distribution n(z) for one particular wavelength is a function of z 
only. Kamb’s expression for the sliding velocity is (Kamb, 1970, p. 703) 


(1+ me2p2\(n-1)/2 


Ub = Aqnt2en—lpnti 


\ATE (3.30) 
or 
_ 2411 + (ee /2)?2)@- DP 2Are 
4enr-1 entl fp’ 
This means that s(0, 0,3) = co(0,3) = 4/e? ~ 0.54, a value that does not fulfill the 
inequality (3.26). Meyssonnier’s results and those of Kamb differ, but not much. In 
light of the assumptions that Kamb had to make the agreement is in fact surprisingly 
good. Kamb’s theory not only gives the correct dependency of the sliding velocity 
on roughness, but the numerical values are also almost correct. Kamb’s non-linear 
theory must therefore be considered to be a successful one.’ 


Ub 





(3.31) 


7There is no doubt that Meyssonnier got the correct result. I have repeated Meyssonnier’s 
calculations based on Johnson's variational theorem (Johnson, 1960), that are easy but involve 
some tedious work, and have got exactly the same numbers. 
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CHAPTER 4 


Characteristics of Flow Close to 
Bedrock Undulations 


In this chapter the flow characteristics of a highly viscous medium flowing without 
friction over a bed consisting of a single sinusoid will be investigated. This is done 
with the help of analytical solutions that are only valid for a linear medium and 
small roughnesses. 


ry 


4.1 Morland’s solutions 


Morland considered the same problem as Nye and Kamb: The creeping flow of 
a Newtonian medium over a wavy bed. He incorporated gravity as the driving 
force of the motion and calculated terms to second order in ¢. For the special case 
of a sinusoidal bed he gave expressions for the pressure field and for the velocity 
components at z = 0, but not for the velocity and the for stress field as functions 
of x and z. Using Morlands results one can calculate without difficulty the velocity 
and the stress field as functions of x and z, although somewhat laborious work is 
involved. Note that only flow in the absence of regelation will be discussed here. 
The corresponding solutions with regelation are given in Appendix A. 


For the velocity field one finds 


ve(z, z) = ut a = (1—2/h)") 


+ 


upak?ze~*? sin kz 
u,(ak)2e~7** (1/4 — kz/2) cos 2kr + O(e*), (4.1) 


+ 


—kz 


v(z,z) = wak(1+kz)e™ coskr 


+ sus(ak)?kze sin 2kx + O(e°). (4.2) 


The sliding velocity is as before 


Tb 


=—. 4.3 
Ub na2k (4.3) 
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The strain rates are given by 


ézz = wUWak zk*e—™ cos kx 
+ us(ak)* k(kz — 1/2)e"*** sin 2kx + O(e%), (4.4) 


én = susk(ak)*(1 — 2/h) 
— wak zk’e7™ sinks 
+ us(ak)? k(kz — 1/2)e~?*? cos 2kx + O(e?). (4.5) 


4.2 Dimensionless form of the flow solutions 


For the following discussion it proves to be of advantage to rescale the various 
dimensional quantities involved and so to recast the equations in a dimensionless 
form. To this end dimensionless vertical and horizontal length scales, denoted by 
capital letters, are defined by 


X :=kz and = Z:= kz, (4.6) 


where the wavenumber k is used as a scaling factor. The velocity field is scaled by 
the sliding velocity 


Vx(X, Z) = uz(z, z)/us and W (X,Y) := v, (a, z) /us. (4.7) 


As dimensionless parameters the slope parameter « and the thinness parameter 6 
will be used. Sometimes Vx will be denoted by U and Vy by V. Eqs. (4.1) and (4.2) 
then become 


2 
Vx(X,Z) = 14+ (1-(1-62)’) 
+ eZe~*sinX 
+ €7e~°4(1/4 — Z/2) cos2X + Ole), (4.8) 


and 
Vz(X,Z) = e(1+Z)e* cosX 
1 
+ 56 Ze sin 2X + O(e°), (4.9) 


where use has been made of Eq. (4.3). 


4.3 Discussion of the second order velocity solution 
with respect to the possibility of extrusion flow 


Does Eq. (4.1) imply extrusion flow, 7.e. does v(x, 2) have a local maximum or 
minimum for z > zo? A necessary criterion for a stationary point of v,(zx,z) is 


Vuz(z,z) = 0. One may start with investigating the variation of Vx(X,Z) with 
respect to X. 
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4.3.1 Variation of Vx(X, Z) with respect to X 


By looking at 


OVx _ > 
aX (X,Z) = eZe~“cosX 
+ €7(1/2— Z)e~*4 sin2X, (4.10) 


one sees that OVx(X, Z)/OX = 0 has as a solution X = 1/2 and X = 37/2 with no 
restrictions on Z.! There is another interesting set of solutions given by 


Ze? 
sin X = e(1— 22) (4.11) 


Since X = 0 and X = 7 with Z = 0 are solutions (these two solution branches will 
now be called (Xo,Zo) and(X,,Z,)) to Eq. (4.11), there are four points where the 
horizontal velocity obtains a minimum or a maximum value at Z = 0 with respect 
to X (and not only two as one might have expected). Inspection of the second 
derivative of Vx(X,Z) with respect to X shows that X = 1/2 and X = 37/2 are 
points of local minima, and that there are local maxima at X = 0 and X = 7a for 
Z = 0. The horizontal velocity does therefore not attain its largest value at the 
highest point of the crest of the sine profile, but at the point of maximum slope. 
This is a second order effect, as can be seen by looking at Eq. (4.8), and is restricted 
to regions in the immediate neighbourhood of the bed where the second order term 
dominates the first order one. There is no first order contribution to the velocity 
field at the bed. In Fig. 4.1 the value of X according to Eq. (4.11) is depicted as a 
function of ¢ and Z. The figure shows that the Z values never get large (Z cannot 
get larger than 1/2), and that as Z increases, X approaches the value X = 1/2 
rather rapidly. At this point a further increase in Z will increase the right-hand side 
of (4.11) to greater than unity.? Only the solution branch corresponding to X = 0 
and Z = 0 is shown. Only in the immediate vicinity of the bed will the maximum 
value of Vx(X,Z) be close to the point X = 0 and X = 7. As an example: for 
€ = 0.1, Z will have to be smaller than 0.03 if X is to be less than 0.1/7. This is due 
to the fact that for Z = 0 the first order contribution to Vx (X, Z) vanishes and it is 
therefore only the second order contribution to Vx(X,Z), which is responsible for 
the maximum at Z = 0. As soon as Z becomes somewhat larger than zero the first 
order contribution starts to dominate and moves the maximum of the horizontal 
velocity with respect to X to the point X = 1/2. There will then be two (and no 
longer four) stationary points: The point X = 2/2, which — as an investigation of 
the second derivative of Vx(X, Z) with respect to X shows — is a point of maximum, 
and the point X = 32/2, which is a point of minimum. Except for this second order 
complication, these are the only X values where one can expect stationary points. 


1Tt is to be understood that because of the periodicity of the sine function an integer multiple 
of 27 can always be added to the values of X although it will not be explicitly so written. 

?As Z increases further the right-hand side of Eq. (4.11) has a local minimum at Z = 1, 
obtaining the value —e/e, and since the absolute value is larger than 1, there are no solutions for 
Z > 1/2. 
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Fig. 4.1: Stationary points at X = 0 and X =7 

Lines of constant X/7 values according to Eq. (4.11), as a function of e and Z. X, Z ande 
are related through X = arcsin{ 725575} /". Each line corresponds to one particular value of 
X/m. The dashed line corresponds to X/m = 1/2. Only above this line is there a solution to 
Eq. (4.11) for some values of Z and €. There are two solution branches to this equation. The 
(Xo, Zo) branch, which includes (X, Z) = (0,0), is shown. The other branch can be obtained 
by the mapping X 1» 1 — X. Notice how rapidly X approaches 7/2 as Z increases, and that, 
for all reasonable values of €, Zp) remains close to the bed. 


4.3.2 Variation of Vx(X, Z) with respect to Z 


Differentiating Vx(X, Z) with respect to Z and setting the resulting expression equal 
to zero gives 


1 
1-6Z= z(4 - le? sin X + (1 — Z)e~?4 cos2X. (4.12) 


The interesting cases to considered are X = 7/2 and X = 3m/2, but notice that 
Z=0, X =0 and X = 7 are also solutions to (4.12). These points, which are 
situated at the bedrock interface, are saddle points. Here the horizontal velocity has 
a maximum with respect to X but a minimum with respect to Z. This minimum 
is situated at Z = 0, which forms the boundary of the medium. The existence of 
these points is, as said before, a second order effect. 


Let us begin with the case X = 1/2 and see if there is a solution to the resulting 
equation 


1 . 7 
1-62 =-(Z- Ie 24 (Z-1)e4, (4.13) 
=:L(n/2,Z,6) 
=:R(1/2,Z,e) 
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a £ = 0.100 


~~ e = 0.138 


R&L 





Fig. 4.2: L(m/2,Z,6) and R(/2, Z,e) as functions of Z. 

The solid lines represent the left-hand side of Eq. (4.13) for a few different values of 6, and 
the dotted lines show the right-hand side of that same equation for different « values. For a 
given 6 value, if € is small enough there will be at least two solutions. If 6 > 0 a third solution 
close to or above the surface (at the surface Z = 1/6) will exist. For every 6, € has to be 
smaller than some critical value if there is to be at least one solution. There are no solutions 
for Z < 1, and as é€ increases the first solution (going along an € curve from left to right) 
moves away from the bed, and the second solution towards the bed. The first solution is, as 
explained in the text, a point of relative maximum of Vx(X, Z) called U5", and the second 
solution is a saddle point called Us%5"'°. 


This is a non-linear equation that does not have a solution in closed form. By 
plotting the left-hand side (L(7/2, Z,6)) and the right-hand side (R(1/2, Z, €)) sep- 
arately, as is done in Fig. 4.2, one sees that there will be a solution to (4.13) if the 
right-hand side, for at least one value of Z, becomes greater than 1 — 6Z. This will 
happen if ¢ is less than some particular value, that will now be called €criticai(7™/2, 4). 
There will only be stationary points at X = 1/2 for some 6 if € < Ecriticai(7/2, 4). 
To determine €-ritical(™/2,6), write € as a function of Z and 6 


Z-1 
(1 — dZ)e2 + (1 — Z)e~2 
Then €critical(7/2,5) can be found by maximising ¢ given by (4.14) subjected to 
Z>Oand0<6 <1. Differentiating «(Z,5) with respect to Z gives 


Ge(Z, 6) _ (2-2 + (2? ~ Z ~ 1)d)e? + (—2? +22 ~ Nen™ (4.15) 


OZ ((6Z — 1)e2 — (Z + 1)e7)? 


By solving de/OZ = 0 numerically, Ecriticai(™/2,6) can be calculated as a function 


e(Z, 6) = (4.14) 





37 


of 6.3 As an example let’s assume 6 = 0, which according to Eq. (4.15) gives 
(Z — 1)?e~?4 = 2 — Z, or — as can be found by doing a few numerical iterations 
with Z = 2 as a starting value — that Z ~ 1.9816906, which gives (c.f. Eq. (4.14)) 
€ & 0.1378839.4 Another interesting limiting case to consider is 6 = 0 ande +0. 
Then, as can be seen by looking at Eq. (4.13), Z must obey (Z — 1)e~? = 0, which 
gives Z = 1 or Z = +00. For infinitely thick “glaciers” with € truly small (let’s say 
€ < 0.05), there will therefore only be one stationary point for X = 7/2, situated at 
Z 1. As « gets larger this stationary point moves upward towards Z ~ 1.9816906 
and another stationary point appears close to the surface which moves downwards 
towards Z ~ 1.9816906 and disappears as e becomes equal to 0.1378839.... At 
this value of €, which will now be called €critical(™/2, 5) (Ecritical(™/2, 0)  0.1378839), 
these stationary points coincide and disappear. Accordingly the point where these 
two stationary points meet as € — Ecritica! Will be called Zeriticai(7/2, 5). We have 
Zeritical(7/2,0) = 1.9816906. An inspection of the determinant of the Hessian matrix 
of Vx(X, Z) at the point X = 7/2 shows that the stationary point below Zeritica(7, 6) 
is a point of relative maximum, so that it will now be called UM. The stationary 
point situated above Zeriticai(7/2,5) is a saddle point; as a function of X, Vx(X, Z) 
has a maximum there, but a minimum if the variation with respect to Z is considered. 
This point will be called Us%9"". 


Ecritical(™/2, 6) is depicted in Fig. 4.3. The figure was obtained by solving 
de(Z, 5)/AZ =0 (4.16) 


numerically for different values of 5, Ecritical(™/2,5), subject to 0?«/0?Z < 0 and 
Z > 1. The figure shows that the value of c, above which the horizontal velocity 
has no stationary points at X = 7/2, gets larger as 6 increases. 


The variation of €critica(7/2,6) and of 6 as a function of Zeiticai(7/2,6) is given 
in Fig. 4.4. As e increases from 0 to €Ecriticai(7/2, 4), m2, goes from Z = 1 to 
Z = ZLeritica(T/2, 5). Zeriticas(7/2,5) increases somewhat with €criticai(7/2,6). For 
every € less than €critical(7/2,6) there will be at least two stationary points. (The 
possibility of a third stationary point will be discussed later.) One of these two 
stationary points is a saddle point and will be called Ussd4e. ya is situated 
above Zcritical(™/2, 6). The other one is a point of relative maximum and will hence 
be called U3‘. Urs" is situated below Zeritical(™/2,5). As € increases U3 moves 
upwards toward Zeritica\(™/2, 5), US moves downwards toward Zeriticai(™/2, ); 
they meet there for a value of € = €critical(™/2,5) and disappear. In Fig. 4.5 the 
position of Us and Us#5"'° as a function of € for two different 5 values is shown. 
It can be seen how these two points move along Z toward each other as € increases 
— Usagile starting from Z = 1/6 and Un starting from Z = 1 — meeting at 


*There are two solutions to 0e/8Z = 0, one corresponding to 07¢/07Z > 0 and another which 
gives 0*«/9?Z < 0. If 5 < 1 then the minimum solution (which exists because of the (1 — 6Z) 
term in Kq. (4.14)) will be situated close to zd = 1. For € lying between these maximum/minimum 
values there are three solutions to (4.13). 

“By looking at Eq. (4.13), setting 6 = 0 and ignoring the second term on the right-hand side 
with respect to the first one (anticipating that e will be small), one finds easily, by calculating the 
maximum of the resulting equation with respect to Z, that « must be larger or equal to e~? = 0.135 
if there is to be a solution; a good approximation of the more correct value of 0.1378839. 
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Fig. 4.3: Ecritical(™/2, 6) as a function of 6. 

For « and 6 values below the line there will be at least two solutions to Eq. (4.13). One of 
these two solutions corresponds to a local maximum of the horizontal velocities above the 
peak of the sinusoidal bed and ts called U7"3*. The other solution is situated above Us 
and corresponds to a saddle point of the horizontal velocities, where v, has a local maximum 
with respect to X but a local minimum with respect to Z, and is called ust. If, for a 
given 6 the roughness of the bed corresponds to a € lying on the curve, there will only be one 
solution to Eq. (4.13). This solution corresponds to the situation were the two points Un 
and Ula coincide and disappear. For e larger than €critica! for some given 5, which in the 
figure corresponds to € and 6 values lying above the line, the horizontal velocity will have no 
local maximum or a saddle point at x = 7/2. Note that €.riticai(7/2, 45) increases as a function 
of 6. There is a limiting value of €criticai(7/2,0) = 0.1378839 as § — 0. 


Z = Zeritica(7/2,5) and € = Ecritica(7/2, 5). It can also be seen that Ecritica (7/2, 6) 
increases as 6 becomes larger (this can be seen better in Fig. 4.3). 


At some point the curvature of R(7/2,Z,¢) will change sign and the solid lines of 
Fig. 4.2 can cross the dotted ones not only twice but three times, giving rise to 
the third stationary point. Calculating d*R(1/2, Z,c)/dz? and setting the resulting 
expression equal to zero gives the position of the inflexion-point of R(7/2, Z,<) 


Z-3 3, 
= —_—__ 4.17 
= aga (4.17) 
so that the inflexion point will be at Z slightly greater than 3. The third stationary 


point, at which the streamwise velocity will be called us fo» must therefore be above 


Z =3s0 that 6 has to be less than 1/3. Use will be closest to the bed if Eq. (4.13) 
is fulfilled where the slope of the left-hand and right-hand sides are the same.” This 


5Then, at the point where the dotted lines of Fig. 4.2 touch the solid ones, the slope of the two 
lines will be the same. 
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Fig. 4.4: Ecriticat(7/2,6) and 6 as a function of Zeriticai(7/2, 4). 
Zeriticai(/2, 6) is the point where the stationary points Um and uss meet and disappear 
as € goes tO Ecritical(™/2,6). For € < Ecriticai(7/2,6) the stationary point nye. is situated 
below Zeriticai(™/2,5) and the stationary point U%934'* above Zeriticai(™/2, 4). critical and 4 
are connected, as is shown in Fig. 4.3. 
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Fig. 4.5: The position of Um and Usage as a function of e for 6 = 0 and 6 = 0.1. 

Zeritical(™/2,5) and Ecritica(7/2,5) are the points where the the lines are vertical. The 
branch above Zoriticai(7/2,5) corresponds to Ue and the one below to Um. For 
€ = Ecriticai(™/2,6) there are no Um and Uss3‘'© points. These solutions were found by 


solving Eq. (4.13) numerically. 
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must happen for Z > 3 if there is to be a third stationary point (Uo) ) at all 
OL(1/2, 2,6) _ OR(1/2, Z, é) 


Un/2, 2,8) = Rn/2, Ze) and = C (4.18) 
Solving these two equations simultaneously for € and 6 yields 
Z-Z-1 
é e (4.19) 


~ 972-22 —-1—-e2° ’ 


and 
pa 0 (2? 22 +1) - 242 (4.20 
- Z2-Z-1 am) 
subject to Z > 3. 


In ne: 4.6 « and 1 — 6Z, where Z in this context stands for the Z coordinate of 
yee jar are plotted as functions of 6, according to Eq. (4.19) and Eq. (4.20). If 6 goes 


to zero, UUP 7/2 B0eS to +00, so that 1 — bU,F, > 0, ie. Up — 1/6, which means 
that ULF /2 18 at the surface. This behavior is expected because as 6 — 0 the glacier 
becomes infinitely thick and the surface reacts only to the average form of the bed, 
which means that close to the surface there is no horizontal velocity variation of the 
vertical velocity component (OVz/0X = 0), and because of the free surface condition 
no shear stress and therefore no shear strain rate, which means that OV; /OZ is zero 
also. The interesting thing about Fig. 4.6 is the fact that this point of no vertical 
variation of the horizontal velocity, which is normally at the surface, can move, for 
some value of 6 and €, somewhat below the surface. 


As ULE » approaches 3 the slope of 1 — 6U, = becomes infinite; if that curve were to 
be followed further it would have a sieeatie slope. This continuation of the curve 
does not correspond to UV? but to Use U, ve does no longer exist. For rather 
large values of 6, larger than about 0.2 es = 0. 20198 corresponds to Z = 3 according 
to Eq. (4.20)), and ¢ larger than 0.35, the vertical variation of Vx will therefore not 
be zero at the surface Z = 1/5 (9Vz/0X # 0) and one must expect the surface to 
start to deform because of bed undulations. This theory can hence not be used for 
6 > 0.2, since the existence of a free surface at Z = 1/5 was assumed. 


The other possible X value for a stationary point of Vx(X, Z), in addition to X = 
m/2,is X = 37/2. Inserting X = 37/2 into Eq. (4.12) gives: 


1-6Z =e-7(1— Z)(1/e —e77) (4.21) 


The left side will always be positive, as will the term 1/e—e~*. For Z = 1 the right- 
hand side is equal to zero, above Z = 1 it is negative and no solution is possible, 
while below Z = 1 it is positive and there will always be a solution if 1 < 1/e — 1, 
or if € < 1/2, which means that €criticai(37/2,5) = 1/2. For ¢ — 0 this solution will 
be situated at Z = 1 and for « = 1/2 it will be situated at Z = 0 no matter what 
the value of 6 is. Again, by looking at the determinant of the Hessian matrix, one 
finds that this stationar¥ point is a point of relative minimum, and it will be called 
Ug). Eq. (4.21) can be solved for e: 


1-Z 


= qc ene (4.22) 
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Fig. 4.6: « and 1—6Z, where Z denotes the normalized vertical position of UV? 7/21 a8 a function 


of 6. Usp is defined as the point of the horizontal velocity above Ga where its vertical 
variation is zero. The coordinate system is scaled in such a way that 6Z = 1 corresponds to 
Un being situated at the surface of the glacier. Since the vertical variation of the horizontal 
velocity is always zero at the surface, which follows from the fact that at the surface the shear 
stresses must be zero, one expects Ue to be found at the surface, and 1 — 6Z to be zero. If 
1—6Z £0, Uns is no longer situated close to the surface, and one must expect the surface 
to start to deform because of the bed undulations. Morland's solutions can then no longer be 
used. Note that the vertical position of Oe. for a given 6, also depends on «. The curve for 
1 — 6Z corresponds to the “worst case”, in the sense that they were calculated for the value 
of e that causes, for a given 6, Uh to be situated as far below the glacier surface as possible. 
This value of € can also be seen and is represented by the « curve. 


As 6 varies the position of Ea changes somewhat. This is depicted in Fig. 4.7, 
which shows the position of Us as a function of ¢ for two different 5 values. A 
value of €critical(37/2,6) = 1/2 is rather large, and it is not clear if this prediction 
can be trusted since the perturbation approach is only valid fore < 1. 


The existence of the stationary points U3", Us9$""* and Uy, shows that there 
will be two regions of extrusion flow. One is at t= = m/2, which extends over the 
region that lies between the saddle point (U%$""*) and the maximum point (Uy"3*), 
and another one at X = 37/2 that eed from bed towards the point of local 
velocity minimum (Ujn"?,). The extrusion flow at X = 7/2 will only be found for 
E < Ecriticat(™/2,5), and the extrusion flow at X = 37/2 will only be found for 


é< 1/2. 


To get an idea of how strong the velocity maximum at X = 7/2, is with respect to 





°For a given € value and 6 ¥ 0 there will be two solutions for Z of Eq. (4.22). Only one of them 
will be situated below 1/6 (which is the maximum value of Z) and it is depicted in Fig. 4.7. The 
other one will always be somewhat above 1/6 and is therefore of no relevance. 
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Fig. 4.7: The position of Uy'7, as a function of « for 6 = 0 and 6 = 0.2. Ursin is the 
local minimum of the horizontal velocities above the trough of the sinusoidal curve (where 
x = 3/2). For e = 0, Un is situated at Z = 1. As € increases it approaches the bed, and 
disappears at ¢ = 1/2. Note that 6 has almost no effect on the position of Us". 


the velocity at (X, Z) = (7/2,0), R™*(e, 6) is defined as 
Une —~ U(n/2,0) 


R™*(c,6) = “Oe Gia) (4.23) 
For 6 = 0 and e > 0 
R™*(¢,0) = eet Vet e(l/e’ +1)/4 (4.24) 


1-7/4 


or R™*(e,0) = e/e (6 = 0, € 40). This equation cannot be used for € values close 
tO Ecritical(™/2,0) because the variation of the position of U7’5* with respect to € has 
been ignored. 


In a similar way R™"(e, 5) is defined as 
Sayy — U(3m/2, 0) 


U (37/2, 0) ee 


R™(e, 6) := 
R™"(e,5) gives an idea of how strong the velocity decrease at X = 37/2 is, with 
respect to the velocity at the bed. 
Kq. (4.8) gives: : 
; —lj/e+e+(l1/e*+1)/4 
min — | 4.26 
R™(e,6) =e 1 e/4 (4.26) 


There is a minimum at: 


e = (Se? +1— V25e4 + Ge? + 1)/e (4.27) 
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Fig. 4.8: Vx as a function of X and Z for e = 0.01 and 6 = 0. The maximum at X/m = 1/2 
and Z = 1 is U3". The minimum at X/n = 3/2 and Z = 1 is Uf2,. For € + 0 the position 
of Uns" and Usr i approaches Z = 1. With increasing c, U3" moves upwards and disappears 
for € = Ecritical: Un /2" on the other hand, moves downwards with increasing € and disappears 
for ¢ = 1/2. 


or € = 0.144. For this value the extrusion flow is the largest. This value is, however, 
only approximately correct since here again the variation of U377, with respect to € 
has been ignored, and its position at Z = 1 for e = 0 was used. 


It is illustrative to look at the velocity field given by (4.8) for some values of ¢ and 
6. In Fig. 4.8 Vx(X, Z) is shown for ¢ = 0.01 and 6 = 0. UMS at X = 1/2 and 

jnyp at X = 31/2 can be seen but not, of course, U375""° since it is at Z = +00 for 
6 =0. Notice that the X is in units of 7. 


The solution of Nye and of Kamb has, as said earlier, regions of extrusion flow. For 
the stationary points this solution gives, for the vertical position of the maximum 
and the minimum points, the value Z = 1. The solution of Nye and of Kamb 
has no saddle point. These findings are reproduced by Morland’s solution in the 
limit ¢« — 0 as is to be expected. For € values as large as €critical(7/2,0) there 
are profound differences. Nye’s and Kamb’s solution makes no distinction between 
the maximum and minimum points regarding their appearance/disappearance with 
respect to different slope parameters values. Morland’s solution predicts instead 
the disappearance of the maximum point, UM, and the saddle point, Usage, at 
€ > Ecritical(7/2,6). The minimum point, Jays remains for € values up to ¢ = 1/2. 
One sees that Morland’s solution contains regions of extrusion flow, but only for a 
limited range of € and 6 values. 


It is important to notice that the influence of the wavy nature of the bed on the 
flow is not largest at the bed but at some distance above it. Typically this distance 
will be about 1 (or z = 1/k). Often the deformation of a bore hole is used to get 
information on the rheological behavior of the ice. One must carefully interpret the 
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Fig. 4.9: Vx as a function of X and Z for « =0.1 and 6 = 0. 
Um has moved upwards and Usn'7, downwards with respect to Fig. 4.8. The point U*754'* can 
also be seen. As € increases further U3 and Us9g‘!* move towards Z = Zeritica(™/2,6) © 
1.9816906, which they reach for € = €citicai(7/2,0) + 0.1378839. Simultaneously Up, 
moves downwards and reaches Z = 0 for € = Ecriticai(37/2,6) = 1/2. 


data from the lowest portion of a bore hole since the exact form of the bedrock, which 
is usually not known, can have a large effect on the flow. Fig. 4.9, for example, shows 
that the perturbed flow (the second and the third terms of Eq. (4.8)) dominates in 
the region Z < 3, the gravity-driven plane flow (the first term of Eq. (4.8)). 


In a numerical study Meyssonnier (1983) observed the appearance and the disap- 
pearance of Um and Us#3"* but did not mention U3". He did not include gravity 
in his model. The driving force of the motion was a given horizontal velocity at. some 
distance above the bed. In this respect his calculations correspond to the problem 
solved by Nye and Kamb. The maximum point disappears because of the super 
position of a linear velocity profile on the Nye and Kamb solution. 


4.4 A physical explanation for extrusion flow 


Upon reflection it becomes evident that the vertical contraction and expansion of 
the ice close to the bed is responsible for the extrusion flow. At some distance, let 
us say at z = 2), sufficiently far above the bed the ice moves parallel to the mean 
bed slope. For kx = 1/2, a high-pressure zone develops above the bed, which causes 
a Poiseuille flow, superimposed on the gravity-driven plane flow (GDPF) solution. 
The maximum of the Poiseuille velocity profile is at (2; — 29)/2, and if its decrease 
above that point is faster than the increase of the GDPF velocity profile, a velocity 
maximum will be found. Since the influence of the bed profile on velocity field is 
(because of the factor e~**) limited to a zone of height proportional to 1/h, 2 will 
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Fig. 4.10: Vx as a function of X and Z for « = 0.1 and 6 = 0.1. 
By comparing this figure with Fig. 4.9 the influence of 6 can be seen. 











Fig. 4.11: Vx as a function of X and Z for ¢ = 0.5 and 6 =0. The points U™3* and Us35i'° 
can no longer be seen, and point U;n'7, is at Z = 0, so there is no extrusion flow. 


be proportional to 1/k and one will expect the position of this maximum to also 
be proportional to 1/k. As a a matter of fact (2.23) has a maximum at z = 1/k 
for ka = 7/2. At kx = 3/2 the ice is expanded vertically and the Poiseuille flow 
profile reverses, causing a velocity minimum, again at z = 1 /k as said before. 


Nye’s and Kamb’s solution ignores GDPF and actually expresses nothing but this 
contraction and expansion due to the wavy nature of the bed. The velocity minimum 
and maximum therefore never disappears no matter how é is varied. If on the other 
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hand the GDPF is present it induces a subtle interplay between the Poiseuille flow 
and the GDPF. GDPF must be involved if extrusion flow is to develop. Therefore 
only certain € values give rise to this interesting flow behavior, as was shown in 
Sec. 4.3. 


4.5 Higher harmonics 


Another interesting feature of Eqs. (4.1) and (4.2) is the doubling of the sine-wave 
frequency. These second harmonics may be somewhat surprising since one might 
expect the reaction of a linear medium to be of the same frequency as the applied 
perturbation, a fact well known from linear system theory. The perturbation is, 
however, in this case not a small term added to the fundamental differential equa- 
tion, as in the usual perturbation theory, but rather the boundary condition that. is 
perturbed and that is a non-linear perturbation. 


A simple argument shows that one is to expect other frequencies than sinkz to 
appear in the solution, unless ¢€ is so small that the first order solutions (cf. Sec. 2.3.1) 
are valid, and that they will become more pronounced as € becomes larger. 


Let us suppose that there were no higher harmonics in the expressions for u and w, 
for all values of ¢. That is 


u=upt+eéssinkx and w=¢,coskz, on z= 2p, (4.28) 


where ¢) and ¢; are some unknown constants. It follows that 


ees a ee (4.29) 
w €,coskr 
Using the exact boundary condition (2.15) and z) = asin kz one obtains 


u 1 

w  akcoskz ae) 
at the base. Comparing Eq. (4.29) and Eq. (4.30) shows that ¢, must be equal to 
upak and that é has to be zero if expression (4.28) is to be true. But this must be 
true for all « values because Eq. (4.30) is always valid. On physical grounds Cg = 0 
can, however, be rejected; for high € values v, will certainly not be a constant on 
z = 2, which means that the starting-point Eq. (4.28), must be incorrect. The 


velocity will therefore not be a single harmonic.’ 


4.6 The stress field 


Using results from Morland (1976a) the stress field can be calculated. My result is 


One = T(1 — z/h) + 2nk* wa {kal sz — 1/2)e~*** cos 2kx — kze~** sin kr} + O(%), 
(4.31) 


"For small corrugation values é) = 0 and é; = usak may become approximately correct. As 
a matter of fact this is what Equations (2.23) and (2.24), that only show the dependence on the 
fundamental harmonic, imply for z = 0. The expressions (4.1) and (4.2) that show the existence 
of the second harmonic do not, of course, conform with these values of ¢p and ¢). 
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Ors — Orr = —4uyk?a { ze™* coskr + ka(kz — 1/2)e~?*? sin 2ka} + O(e?), (4.32) 


The expression for the second invariant of the deviatoric stress tensor is particularly 
interesting: 


?o= P(1—2/h)y’ 


- 4rp=(1 — z/h)e~*™ sin kx 
+ arb{(=)Pe™ + (1 — z/h)(kz — 1/2)e7?** cos 2k} 
Bry we (kz — 1/2)sin kr 
+ 4Arpe***(kz — 1/2)? + O(e°). (4.33) 


This equation should be compared with the corresponding Eq. (2.28) of Nye and 
Kamb. In contrast to (2.28) it shows that 7 is dependent on zx. It can readily be 
shown that if z/h < 1, Eq. (4.33) has stationary points at (kr, kz) = (pi/2, 1), which 
is a saddle point, and (kz, kz) = (32/2,1), which is a point of absolute maximum. 
These stationary points do not disappear at certain corrugation values as did U7", 
Usass¢ and Us’). Notice that r, for all values of z, increases from the bed upwards, 
attaining its largest value at kz = 1. To get an idea of how large the variation of r 
with respect to x is, consider the ratio: 


__ tlkr = 1/2,kz =1) 

“(kz = 30/2, kz = 1) 
which is the ratio of the maximum and minimum values of 7 as a function of z. 
Putting (4.33) into (4.34), and ignoring z/h with respect to other terms, gives 


_ [1 —2/e? — 1/e* + 4(1 — 1/e)/(ee) 
a 4(1 — 1/e?)/(ee)1 + 2/e? — 1/e4 e33) 
0.71lée + 1.272 
RS Vo 7aBe 1.072 eae) 


For « = 0 there is of course no variation of t with z, 7.e. R = 1. A corrugation 
value of ¢ = 0.1 gives R ~ 1.12 or 12%, and € = 0.3 gives a variation of about 
about 42% in 7 as a function of xz, for kz = 1. This shows, without any doubt, that 


Kamb’s main assumption in his development of a non-linear theory of sliding, €1; 
being independent of z, is only approximately valid. 


(4.34) 


or 


Kamb never claimed that €;; was entirely independent of x but merely that the flow 
is dominated by the z dependence of 7. For small enough ¢€ values this will be true. 
The asymptotic « dependence of uy for « > 0 is given correctly by Kamb’s theory, 
although the numerical values are not exact, as has been shown on page 32. 


All of the above equations for the velocity, the stress field, etc. are only valid as long 
as no bed separation takes place, that is as long as the pressure at the bed remains 
positive. The pressure field is given by 


tan a . : 
P(x, z) = pa + prghcosa {1 —2z/h+ = ee" cos kz + kae~?** sin 2kx)} + O(e°) 


(4.37) 
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where a is the mean bed slope and p, the atmospheric pressure. 


The position of the points of maximum and minimum pressure along the bed are 
found by solving Op(z, z)/0z = 0 with z = 0, which gives 


= sin kx 
~ cos 2kx’ 





(4.38) 


Using a Taylor expansion, one finds that the pressure obtains its minimum value at 
kx = 7 — € + O(e3) and its maximum value at kr = € + O(e3). This shows that as 
e gets larger these points move from the inflexion points of the sine curve towards 


kr = 17/2. 
Requiring p > 0 gives 


E Pa 3 
t <=¢1+ —— O 4. 
ana < | +e (e"), (4.39) 


or to a good approximation: tana < ¢/2, if the atmospheric pressure can be ignored 
with respect to the mean ice overburden pressure. By inverting (4.39) one gets . 





2 
1 a a 
€ > a Seer | (yi eee — (4tana)? 
4tana P1g COS a@ P1g Cosa 
PaSPig°S% 9 tana + O((4tana)’), (4.40) 


which gives a lower limit for ¢, for a given mean bed slope and ice thickness.2 There 
are therefore two restrictions on €. First of all, it has to be small if the results of 
the perturbation analysis are to be applicable, and secondly, it may not be so small 
that a bed separation occurs. The predictions of this theory can only be trusted for 
values of € up to « = 0.5, which means that a has to be less than 14°. Applications 
to ice falls are therefore questionable. 


4.7 Summary 


At this point I would like to summarise the most important results obtained so far. 


The sliding law has been given in a non-dimensional form and it has been argued 
that only ¢, 6 and n will enter it in a nontrivial way. The fact that u, increases 
linearly with \ and a, for a given € independent of n and 6, as shown in Eq. (3.6), 
follows directly from dimensional arguments. 


For asymptotically small ¢ and 6 values, the sliding law is given by (3.23). A 
simple argument was given, which does not prove this result, but makes it physically 
reasonable. There are other theoretical arguments that give the same result (Kamb, 
1970; Fowler, 1979). 


8The solution with a plus sign before the square root in (4.40) can be rejected by noticing that 
it does not give the correct. value in the limit of a — 0. 
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Using the results of Morland (1976a), the expressions (4.1), (4.2), (4.32), (4.31) and 
(4.33) that give the stress and velocity fields as functions of + and z have been 
derived. It turns out that of two of the particularly interesting features of Nye’s 
and Kamb’s solution — occurrence of extrusion flow and the second invariant of the 
strain rates being independent of x — the first is carried on to these more precise 
solutions but the second is not. It must be stressed that this extrusion flow is of local 
character, and should possibly be called local extrusion flow in order to distinguish 
it from extrusion flow encountered earlier in the glaciological literature, which was 
of global nature — extending, at some distance below the surface, through an entire 
glacier — and causing an increase of ice flux. The local extrusion flow does not add 
to the ice flux. 


Criteria for extrusion flow have been given. Close to the bed the horizontal velocity 
field can have at most three different stationary points. If 0 < e < 1/2 there is a 
minimum point (U jaya) situated at X = 37/2 and Z < 1, which moves downwards 
as € gets larger, reaching the bed for « = 1/2. For X = 7/2 there is a point of 
maximum velocity (U5") at 1 < Z < Zoriticat(7/2, 5) and a saddle point (Us35"") at, 
Z > Zeritical(7/2,6) as long as € < Ecritical(7/2, 5). With € increasing the maximum 
point, U3", moves upwards and the saddle point, Us3$""°, downwards until they 


meet at Z = Zeriticai(7/2, 6) for € = Ecritical(4/2, 6) and vanish. 


A simple argument was given, which shows that the velocity will in general not be 
a simple harmonic at z = zo, except for truly small values of €; a feature that may 
at first be somewhat surprising but is evident in expressions (4.1) and (4.2). At 
large corrugation values this frequency doubling may be expected to become more 
pronounced. 


CHAPTER 5 


Testing of the Correctness of 
Numerical Results Obtained with the 
| FE Program MARC 


The general purpose FE Program MARC (MARC, 1992) was used for all numerical 
calculations of flow velocities and stresses presented in this work, including calcula- 
tion of flow over a sinusoidal bed. In this chapter various tests done to estimate the 
errors involved with the numerical calculations are described, and error estimates 
given. 


The results obtained by the FE code were tested by comparing them with analytical] 
solutions for some analytically solvable problems. For the numerical calculation 
of flow over a sinusoidal bed an additional type of testing was done by comparing 
results of runs with different values of A, a, h and a but same « and 6 values together. 
This latter type of testing gives an estimate of the relative errors that are caused by 
differences in grid density. Comparison with analytical results, on the other hand, 
gives an estimate of the absolute errors involved. It is important that both types 
of testing are done, since the problems that can be solved analytically often possess 
some kinds of symmetries that are not representative of real glaciers; comparing 
numerical results with them does therefore not necessarily test the programs in a 
general way. The discretization error was also estimated by solving the same problem 
using several meshes with different numbers of elements. 


5.1 Comparison with analytical solutions 


When comparing calculated results with analytical ones, one should not use a prop- 
erty of the analytical solutions as a boundary condition (BC) for the model. For 
example, calculating the flow of a gravity-driven plane-flow using the analytical so- 
lution for the velocity profile or the stress variation along the edge as a boundary 
condition for the model will give estimates of the numerical errors which are too low. 
A better approach would be to use periodic boundary conditions. Results obtained 
with MARC were compared with analytical solutions for five different. problems. 


5.1.1 Gravity driven plane flow down an inclined plane 


The effect. of using velocity or pressure as a boundary condition, and the use of 
periodic boundary conditions on the accuracy and computer time required, were in- 
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vestigated. With 8 linear elements in the vertical direction the maximal discrepancy 
between numerically and analytically calculated velocity for n = 3 was 0.04% if the 
velocity was given as a BC, 0.06% if the pressure was given as BC, and 0.1% for 
periodic boundary conditions. For non-linear rheological behavior the effective vis- 
cosity distribution was repetitively calculated until the maximum change in velocity 
was less than 0.01%. The iteration procedure converges much faster if the velocity is 
given as a BC, thus confining the velocity variation, than if the pressure or periodic 
boundary conditions are used. 


5.1.2 Closure rate of a circular hole and the flow between two rotating 
coaxial cylinders (Couette Flow) 


The results concerning the numerical accuracy of the calculation of the closure rate 
of a circular hole, are similar to those obtained by studying gravity-driven plane flow 
down an inclined plane. It turned out that for a given mesh the numerical results 
for flow between two rotating coaxial cylinders were more exact (about an order of 
magnitude) than those for the closure rate of a circular hole with some pressure p; 
within the hole and some pressure p2 at the outer boundary. This shows that the 
accuracy of the numerical results obtained with the same mesh can vary greatly for 
two seemingly similar problems. 


5.1.3 Gravity driven flow down a channel of circular cross-section 


The velocity profile for a semi-circular channel of radius R of uniform cross-section 
and uniform slope a is 


oemen(3) AO-(@"). 


where 7 = p;gRsina, and r = \/y? + z2 (Nye, 1965; Hutter, 1983). The discharge 
q is 





2 
l| 


a dé ie dr u(r, 6) 


w/2 


a 3) (3) ne 4 





If the velocity is scaled with 2AR7{' and the discharge with 2A7? R? the dimensionless 
velocity Up and discharge Q are given by 


1 
U) = ———~ and Z 


2(n +1) es 2°41 (n + 3) ee) 


These scalings are the same as those used by Nye (1965).! 





'The A that Nye uses is twice the A used here, i.e. Anye 1965 = 2A. 








= 1.0 
2 08 
+ 06 re 
UO — 
R04 c 
= tJ 
E 0.2 
= 0.0 wT 9 nn PR eer 
-100 -—50 0 50 100 -100 -50 0 50 100 
Distance from center Distance from center 
Fig. 5.1-a: Normalized velocity profiles Fig. 5.1-b: Percentage difference between 
along the z-axis and the y-axis for n = 1. analytical and calculated v, values for n = 
Detailed information is given in the text. 1. The long dashes represent the +0.15% 
range within which most of the values lie 
within. 
- 1.0 
© 0.8 
Z 0.6 es 
N 0.4 c 
Bb uJ 
E02 
= 0.0 Pusey soap eUsEEssyuerte = 
-100 -50 0 5 100 -—100 -50 0 50 100 
Distance from center Distance from center 
Fig. 5.1-c: Normalized velocity profiles along Fig. 5.1-d: Percentage difference between 
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Figs. a and c show the velocity variation of all the three spatial velocity components along a 
profile that runs from a point at the bottom of the center of the channel (y = 0 and z = —R) 
upwards along the z-axis towards the centerline (y = 0 and z = 0), and then sidewards across 
the surface along the y-axis toward the margin (y = R and z = 0). Calculated values are 
represented by symbols. v, is shown as a plus sign, v, as an asterisk, and v, as a diamond. 
The solid line that goes through the plus signs in Figs. a and c is the analytical result for v_. 
The perceptual difference between analytically and numerically calculated v, values is given in 
Figs. b and d. All the figures refer to the 10 x 10 discretization. Since the transversal velocity 
component v, was fixed (= 0) along the longitudinal central plane of the channel, v, is not 
shown in the left parts of Figs. a and c. 


For the calculations two FE discretizations where used, one with 10 x 10 elements 
in the cross-section”, and another finer discretization with 20 x 20 elements in the 
cross-section. The meshes were generated by an automatic 2D mesh generator. An 
example of the meshes created in this way can be seen in Fig. 5.2. The domain of 
investigation was three-dimensional, but was reduced to two dimensions by imposing 


210 elements subdivision in each spatial dimension, giving a total of 100 elements. 
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Uo Q 
theoretical value 0.03125 0.032725 
10x10 mesh 0.03138 (0.42%) 0.03321 (1.5%) 
20 x 20 mesh 0.03126 (0.03%) 0.032703 (0.06%) 


Table 5.1: Semi-circular channel. Comparison of analytical and numerical values for Up and Q 
for a semi-circular channel with n = 3. Values in the brackets stand for the percentage error 
in Up and Q. 


periodic boundary conditions for the velocities. No assumptions about the stresses 
or the velocities at the “cut off” boundaries were made. Only one half of the semi- 
circular channel was modelled. 


The difference between analytical and numerical values can be seen in Table 5.1. 
‘The error is considerably smaller for the 20 x 20 discretization, but even for the 
10 x 10 discretization quite acceptable. It can be concluded that a mesh with about 
10 elements over the depth of a glacier gives sufficient precision for most practical 
purposes. Since a discretization with more than 10 elements over the depth can be 
computationally expensive the errors of the 10 x 10 discretization are analysed more 
closely in Figs 5.1-a to 5.1-d. 


Sometimes the flow transverse to the main flow direction is of interest. The FE 
formulation used allows the elements to be exactly incompressible in principle, but 
due to numerical errors this is not exactly so in practice. This causes some transverse 
flow even for infinite channels where no transverse flow should occur. It is important 
to estimate the magnitude of this “flow”, which is linearly proportional to the main 
flow, and which expresses nothing but numerical errors. In Fig. 5.2 the horizontal 
component of the transverse flow (v,), and in Fig. 5.3 the vertical component (v,), 
are shown for a semi-circular channel with radius 100m. The maximum horizontal 
velocity v, was equal to one. This “numerical” transverse flow is surprisingly large, 
up to 3.5 % of the maximum horizontal velocity for v, and up to 2 % for vy. In 
general the errors are, however, much smaller, but this shows how careful one must 
be in interpretating calculated transverse flow results. 


In Table 5.2 computed results for a semi-circular channel for n = 1 to n = 5, where 
n is the power law exponent in Glen’s flow law, and the corresponding errors are 
listed. The errors of the velocity on the centerline do not seem to increase with n 
and are of the order of 0.2 %. The mesh was 3D and the results given are for one 
particular cross-section. Inspection of the results showed the centerline velocity to 
oscillate slightly along the z-direction. Taking this small oscillation into account 
the correct error estimate is 0.22 %. Errors smaller than this are accidental. The 
discharge errors increase somewhat with increasing n. It is not clear whether this is 
a result of the flow field errors or the errors involved with integrating the flow field 
to get the discharge. Since only the resulting error is of interest this question is of 


no concern. 
In conclusion it can be said that the program MARC can be used for 3D non-linear 
calculations of glacier flow where thermal effects are not important. For about 10 
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|| 5.45e-03 
2.32e-03 
-8.15e-04 


-3.95e-03 





Fig. 5.2: Normalized transverse horizontal flow (v,/v,(R = 0)) in a semi-circular channel. 
The figure depicts nothing but numerical errors. If the calculation were exact, there would be 
no transverse flow. The value of A was chosen such that v, along the centerline was equal to 
one. 






0.051980 










Table 5.2: Comparison of numerical and analytical values for the velocity at the centerline and 
the discharge of a semi-circular channel. The numerical values were obtained with the 20 x 20 


mesh. 
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Fig. 5.3: Normalized vertical flow (v./v,(R = 0)) in a semi-circular channel. 


elements over the depth the error is generally less then 0.5% but somewhat larger 
at the edges, or up to about 4%. There seems to be no reason to use more that 8 
to 10 elements over the depth for most types of calculations. 


5.1.4 Gravity driven flow down a channel of parabolic cross-section 


Numerical solutions exist for glacier flow in channels of various cross-sections (Nye, 
1965). Flows through a few parabolic channels having different aspect ratios were 
calculated. The aspect ratio (W) is defined as the ratio of the channel half width 
(bp) to its depth (a,), i.e. W := b,/ay. The suffix p refers to parabolic. 


A program was written that transforms a semi-circular channel into a parabolic one. 
The idea behind the transformation formula used is explained in Fig. 5.4. The node 
p of the semi-circular FE mesh with coordinates (a, y,z), is transformed into node 






(yz) 


Fig. 5.4: Transformation of a semi-circular channel into a parabolic one. The point r with the 
coordinates (y, z) is transformed onto the point r with the coordinates (y,z'). S(y’) is the 
distance from the lower left corner of the parabolic channel (a,) along the lower boundary to 
the point with the coordinates (y ,z ). (y',z_) is defined as the point of intersection of a line 
going through the channels center point (the origin of the coordinate system) and the point 
r with the lower boundary of the channel. The idea behind the transformation i is that r/R 
should be equal to r/R’, and that 6/(/2) should be equal to S(y")/S(b,), where S(b,) is 
the distance from a, along the bottom to 6,. 


p' having the coordinates (z, y’, 2’) of the parabolic mesh by requiring that 





O: ., s(y’) 
7/2 ~ s(t) _ 
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where s(y”) is a function giving the length of the parabolic curve from (y, z) = (0, 2) 
to (y,z) = (y”, z”). x points in the down-flow direction (making a vertical angle a to 
the horizontal), y is horizontal and perpendicular to the flow, and z points upwards. 
Let zo denote the z coordinate of the lower boundary of the parabolic channel: 


zy = ap((y/b,)? — 1). (5.6) 


The function s(y) is given by 


y dzo oe 
[ yia) dj 
zs of Ve Pan Fe ai 
ve 
= ve {wPap +P + 5 — arcsin b,y/,/2a, \ (5.7) 


V2b 
Eq. (5.4) has no closed solution for y". It was therefore solved numerically by finding 
the zero of the function f(y) = 20/m — s(y)/s(b,). This procedure gave meshes that 
were only slightly distorted. An example of a mesh generated in this way is given 


in Fig. 5.5. 
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Fig. 5.5: Parabolic mesh with aspect ratio 2. 


Results of the FE calculations are given in Table 5.3. They are compared with earlier 
results from Nye (1965) wherever possible. The form factor f is defined so that fz is 
the best linear approximation to the shear stress o,, giving the true surface velocity, 
if fz is integrated. A simple calculation shows f to be given by f = ((n+1)Up)/"3 
Nye (1965) gave estimates of the error in Up and the values in Table 5.3 agree within 
the range of errors. 


Calculating the flow of the parabolic channels involved exactly the same steps as 
needed for the semi-circular channels, except for the additional use of the program 
transforming a semi-circular channel into a parabolic one. For error estimates the 
corresponding values in Table 5.2 can be used, as the errors are expected to depend 
mainly on n and to lesser extent on W. 


f and Q as functions of W are plotted in Figs. 5.6 and 5.7 respectively. 


For W — 0 one finds that U = W"+!/(n + 1), and that f = W'"+1)/", and since 
ndf/dW = (n+ 1)W'/" the slope of f(W) at W = 0, for any fixed value of n, is 
zero. On the other hand, if n + oo and W — 0 then the form factor is proportional 
to W, and the slope equal to one. This can be seen in Fig. 5.6 and explains the 
somewhat different slope of the curve for n = 1 with respect to the other curves 


3In general the form factor will depend on the variation of basal sliding across the width of 
the channel. The form factors in Table 5.3 are valid for the special case of no, or uniform sliding. 
Highly non-uniform sliding across a glacier cross-section has been observed by Raymond (1971) 
at the Athabasca Glacier, where basal sliding velocities up to about 40 m/a and at the same time 
negligible marginal sliding was found. Reynaud (1973) showed that this kind of lateral sliding 
variation can for example result from a friction at the bed being proportional to effective normal 
pressure. Further numerical work, that resulted in better agreement with data from Athabasca 
Glacier, has been done by Harbor (1992), who used a sliding law dependent on basal shear stress 


and effective pressure, with the additional assumption of increased friction towards the margin of 
the glacier. 
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Table 5.3: Up is the velocity at the centerline normalized by 2a,Ar7;,", and Q is the discharge 
normalized by 2Aa37;". The form factor was calculated according to f = ((n+1)Uo)'/". The 
numbers in the brackets are from Nye (1965). Calculated values are estimated to deviate less 
than 0.22% from exact ones. 


that at first sight might seem to indicate an error. Notice that the curves in that 
figure result from simple interpolation and that the slope at W = 0 is incorrect. As 
n gets larger, increasingly smaller W values will be needed to make the correct slope 


become evident. 


An interesting feature of Fig. 5.6 is the cross over of the f (W) curves for W . 1.2; 
For larger W values increasing n and keeping W fixed has the effect of decreasing f. 
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Fig. 5.6: Form factor f for parabolic channels as a function of W. Letting W — 0 for a fixed 
value of n gives a slope of zero at W = 0. If, on the other hand, n — oo and then W - 0, the 
slope at W is equal to unity. Note the cross over of the curves. This happens at a value of W 
for which the drag from the bottom starts to be more important than the drag from the sides. 
Since the limiting value of f = 1 for W — oo is approached more slowly for larger values of n, 
increasing 7 for a wide channels has the effect of making the glacier more aware of its finite 
width. For a narrow channel the situation is the same in the sense that increasing n makes the 
glacier more aware of its finite depth. As explained in the text, this can be understood to be 
a consequence of the stress redistribution caused by increasing the value of n. The effective 
stresses then become more uniformly distributed over the bed, and all sections of the bed tend 
to contribute more equally to the overall drag. Increasing n has thus the effect of making a 
valley glacier more aware of the shape of its cross-section. 








Fig. 5.7: Non-dimensional discharge, Q for parabolic channels as a function of W. Q is related 
to the dimensional discharge q through q = 2ATfasQ. 


For smaller W values it is the other way round; f increases with increasing n. This 
shows that for a wide glacier (wide in the sense of W > 1) the effects of the sides are 
felt more closely with increasing n, since the asymptotic value of 1, as W approaches 
infinity, is reached more slowly. Conversely, for a narrow glacier (W < 1) the drag 
from the bottom has a larger effect. on the centerline velocity as n increases. One 
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can therefore say that increasing n has the effect of making a valley glacier more 
aware of the shape of the cross-section. 


The cross over is where the drag from the side starts to dominate the drag from the 
bottom. For an elliptical channel the cross over would hence be expected to be at 
W = 1 but for a parabolic channel this is not necessary true, and as Fig. 5.6 shows 
the cross over is at W > 1. This is of course caused by the curvature of the channel. 
Since the point of the bed closest to the centerline is always the point of largest 
shear stress (Nye, 1965), the location of this point along the transverse section as a 
function of W is important for the relative contribution of the sides to the overall 
drag. For a parabolic channel one can easily show that the point of the bed lying 
closest to the centerline is on the “side” of the parabolic profile, in the sense that a 
straight line from the center to this point then makes an angle that is greater than 
n/4 to the vertical, if W > V2.4 Thus, because of the curvature of the parabolic 
curve, the cross-over point would be expected to lie between 1 and V2, as it does. 


With increasing n the effective stresses across the bed are redistributed and the 
stress distribution becomes more uniform. For a wide glacier this means that the 
relative contribution of the sides to the drag becomes larger. Hence the sides are 
“felt” more strongly and the value of form factor decreases. For a narrow glacier 
the situation reverses. This explains the cross over in Fig. 5.6. 


For W — oo the parabolic channel becomes an infinitely wide channel of uniform 
depth and f + 1 and U > 1/(n +1). 


Nye gives, for the asymptotic slope of the Q(W) curve, 


_ Q_2 


where I(r) = fo’? cos” d0, or 





km 2 - 2 PU/2P(n +3) 


WoW n+2 I(n+7/2) oe) 


where I is the gamma function. This expression is valid for all values of n (and not 
just for integers as is the case for the expression given by Nye (1965)). Hence, for 
n = [(1,2,3,4,5] the asymptotic slope of the Q(W) curve is [0.3048, 0.2032, 0.1478, 
0.1137, 0.09093] respectively. This is in good agreement with Fig. 5.7. The slope 
between W = 3 and W = 4 is for example equal to 0.2173 for n = 2 , and 0.1470 
for n = 3. 

In the limit W —> 0, it can be shown (Nye, 1965) that (n + 2)(n+ 4)Q 7 4W"*?, 
which explains the flattening of the curves in the figure. 


5.1.5 Flow over a sinusoidal bed 


The flow over a frictionless sinusoidal bed was calculated and the results compared 
with the analytical solution given in Sec. 4.1. 


4To see this minimize R? := y2 + 22 with respect to yo, giving y = bp/1— W?/2 and z = 
a,((1 — W?/2) — 1)”, and then solve yo/zo = 1. 
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Fig. 5.8: Calculated horizontal velocity of flow over sinusoidal bed for n = 1. The vertical 
position of the bed zp is given by: z = asinkx. The velocities were normalized by the sliding 
velocity u,. X and Z are non-dimensional coordinates defined as X = kx and Z = kz. The 
peak of the sinusoidal bed is at X/mx = 1/2 and the trough at X/m = 37/2. The values 
of the slope parameter « = ak, and the thinness parameter 6 = 1/(hk), were ¢ = 0.37 and 
6 = 0.05/7. 


In Fig. 5.8 the calculated horizontal velocity field for a/X = 0.015, A/h = 0.1, sina = 

0.01 and pyg = 8.99577 x 10*kgm~?m~? is shown. The domain had dimensions of 

200 times 200 meters, and the FE mesh had 5395 nodes and 5038 elements. For 

the region in the figure the nodal values were used for creating an evenly-spaced 

grid with 200 times 200 points, which was then used for the contouring. The points 
jay and Ur can be seen. 


The velocity field was normalized by the calculated sliding velocity u, = 21.24472 
m/a, which gives s(37/10, 1/(207), 1) = 0.9866, a value close to the theoretical one 
of s(0,0,1) = 1. s(e,6,n) is the sliding function, defined as 

etiku, 


s(e,6,n) := We 
b 





(5.10) 


The mean of the vertical velocity (< v, >) along the bed should be equal to zero, 
but due to numerical errors will never be exactly so. The calculated value was 
<u, >= —1.417 x 10~, or five orders of magnitude less than the sliding velocity. 


In Fig. 5.9 the horizontal velocity field as given by Eq. (4.1) is depicted. Here again 
the velocity field was normalized by the sliding velocity, but the value for u, used 
was not the numerical one, but that according to (4.3). It can be seen that there 
is almost no difference between these two figures by looking at Fig. 5.10, where the 
difference is shown in percentage terms. The difference is small, always within the 
range of —0.06 to 0.12%. It is therefore concluded that the errors involved with 
the numerical solution are less than 0.1%. This error estimate can only be used for 
the calculations of linear viscous media. For non-linear rheology the corresponding 
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Fig. 5.9: Analytical horizontal velocity for n = 1, € = 0.3m and 6 = 0.05/7. Horizontal 
velocity v,(z, z)/us as given by Eq. (4.1). A comparison with Fig. 5.8 shows how similar the 
results obtained by FE calculations are. 


non-linear equation system has to be solved iteratively until a prescribed accuracy is 
obtained. In principle a much better accuracy can be obtained but that may require 
an unreasonably large amount of computer time. 


The iterations were continued until the maximum change in velocity over the whole 
model divided by the maximum velocity was less than some prescribed value, which 
was chosen to be equal to 0.0001 or 0.01 % for all of the calculations. Although this 
is a truly small value the accuracy of the non-linear results will be somewhat less 
than the linear ones. How much larger the errors are is difficult to estimate exactly 
since no analytical solutions are available, but assuming that the only additional 
source of error is the iteration procedure the errors should be less than 0.2%. There 
is, however, another source of error; the effective viscosity distribution can only be 
as accurately represented as the FE discretization allows, with every integration 
point assigned one effective viscosity value. This discretization error does not play 
a role if a linear viscous medium is modelled, since then the viscosity is constant 
over the whole domain. In the non-linear case it must be estimated by solving the 
same problem with different FE meshes. Several calculations using a mesh of 4302 
elements and 4636 nodes were therefore done. The results, compared with the finer 
discretization, differed by less than 0.3%. 


5.2 Discretization errors 


Discretization errors were estiinated by: a) generating several meshes with different 
numbers of nodes and elements, and b) by comparing the nondimensionalized results 
obtained by different meshes characterized by the same nondimensional numbers ¢, 
6 and n, but different dimensional numbers such as a, A, ™ and h. 
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Fig. 5.10: Percentage differences between a numerical and the analytical solutions of highly 
viscous flow over a sinusoidal bed. 

Differences between analytical and numerical results for the horizontal velocity in percents, 
ie. 100(u2(x, z) — vEE(2x, z))/v2(z, z), where a stands for analytical and FE for numerical 
results obtained by an FE calculation, for « = 0.37 and 6 = 0.05/7. One sees therefore the 
percentage difference between Fig. 5.8 and Fig. 5.9. The error is never larger than 0.12%, and 
usually much less. There is a small systematic tendency of the error to be positive for X > 7 
and negative for X < 7. By generating similar figures for different values of the thinness 
parameter 4, it was discovered that this discrepancy is due to the assumption in the derivation 
of the analytical solution that 6 is approximately zero. The numerical solution is therefore even 
more precise than the figure indicates. The analytical solution is defined for z > 0 whereas 
the numerical one is defined for z > z ). No attempt was made to account for this difference, 
which is the reason for the “black” area close to (X,Z) = (1/2,0). 


The lines in Table 5.4 show the percentage difference in s(n,¢,6) obtained from two 
calculations done with different values of the dimensional parameters 4, a, and h, 
which give the same dimensionless ratios « and 6 (r = €/2m and ¢ = 276), and 
should therefore give exactly the same results. The error is almost always less than 
0.2%. For n #1 there is a tendency of the error to get larger as r increases. This is 
most probably due to the fact that as r becomes larger there is an increased stress 
concentration at the peaks of the bed protubances, leading to large variations in the 


effective viscosity over a limited area, which can only be approximately represented 
by the FE mesh. 
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Table 5.4: Comparison of results obtained for different FE grids but for the same dimensionless 
numbers n, r and ¢. Every line of the table shows the values of n, 7, and ¢ for two calculations 
with different values of , a, and h (not tabulated). The last column shows the difference 
in percent of the calculated values for the function s(n,<,6) between these two calculations, 
which should theoretically be zero. 
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CHAPTER 6 


Numerical Calculations of Flow Over 
a Sinusoidal Bed 


Up to this point the emphasis has been on theoretical work done on perfect sliding 
over a sinusoidal bed although some numerical work has also been mentioned. In this 
chapter the flow of a viscous material over a sinusoidal bed for free-slip boundary 
conditions will be analysed. The analytical results that are only valid for slight 
roughness will be extended to extreme roughness and the effect of a non-linear flow 
law on the sliding velocity and on the flow characteristics will be examined. 


6.1 Objectives 


The flow of linear viscous material over a sinusoidal bed for ¢ < 1 and 6 & 0 is 
rather well understood. Eqs. (4.1), (4.2), (4.32), (4.31) and (4.33) gave some new 
insight. In general however, even for the linear case, the dependence of the sliding 
velocity on € and 6 is not known. For a power-law fluid the relationship is still less 
clear since the asymptotic behavior of u, for « — 0 is, as said before, a subject of 
debate. 


Neither is the effect of 6 on the sliding velocity clear. Claims have been made that 
up, depends linearly on 6 (Fowler, 1979) independent of n, but this has also been 
questioned since physically one might expect wp, in the limit of small 6 values, only 
to be dependent on h through its dependence on rT. 


The aim of the numerical modeling was to get answers to the following questions: 


1. How does U, depend on «¢ in the limit ¢ > 0 for 6 < 1? Is Eq. (3.23) correct? 
2. How does us, vary as a function of n. 


3. How does uz, depend on 6? 


4. How exact are the predictions made for the appearance/disappearance of ex- 
trusion flow at certain € values? 


9. Is extrusion flow for n > 1 possible, and how does n effect its magnitude and 
its dependency on ¢€? 
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6.2 Solution method 


The differential equations that are to be solved are the conservation Eqs. (2.4), (2.5) 
together with Glen’s flow law (Eq. (2.1)) given in Sec. 2.1.1. 


The solution procedure applied was the method of finite elements. Commercially 
available FE software, the general purpose program MARC, was used for the calcu- 
lations. This program proved to be well suited, especially because the user can add 
user-defined subprograms to the main code and recompile it. In this manner the 
program can be extended and modified in order to fit the needs of the particular 
person using it. MARC is actually not designed for fluid analysis and Glen’s flow law 
is not implemented. The implementation of Glen’s flow law, after having contacted 
the program developers, nevertheless turned out to be rather easy. One of the strong 
points of MARC is that it allows the calculation of a truly incompressible medium 
and offers iteration algorithms for the handling of systems of non-linear equations. 


A comprehensive testing of the correctness of the numerical results obtained with 
MARC was done. A description of the testing procedure and the results obtained 
is given in Chapter 5. In addition to giving information on the numerical errors 
involved the testing procedure also gave results that are possibly of general interest. 
An example of this is a calculation of form factors for parabolic channels for a 
considerably larger range of parameters than was originally done by Nye (1965). 


The FE element chosen was a four-noded isoparametric element with bilinear in- 
terpolation. Linear elements tend to give better results for non-linear calculation 
although a somewhat larger number of subdivisions may be needed. To enforce 
incompressibility two methods were used: a) selective integration, where the four 
Gaussian points are used for the deviatoric contribution, and the centroid for the 
dilatation contribution, a method also known as the constant dilatation method, 
and b) the Herman variational principle using (v,p)-formulation (mixed formula- 
tion) where a fifth-node for the pressure degree of freedom is added to the center 
of the element. The pressure is constant within the element (it can be shown that 
the pressure interpolation function must be of one degree lower than the velocity 
functions if the resulting equation system is to be non-singular). 


To improve the bending characteristics of the element the assumed strain formula- 
tion was used, where the interpolation functions are expressed with respect to the 
distorted element. 


Due to the non-linearity of the flow law the solution must be sought in an iterative 
way. The method of direct substitution was used, t.e. the old effective viscosity 
distribution was repeatedly substituted by a new one based on the calculated strain 
rates using the old distribution. The iteration was continued until the ratio of the 
maximum change of the velocity at some node in the FE mesh during the last 
iteration, divided by the maximum node-velocity, was less than a specific value, set 
to 0.0001. Experience with other FE programs such as RHEO-STAUB has shown 
that it is sometimes necessary to damp the iteration in order to ensure convergence, 
and to use as a new effective viscosity some weighted average of the old viscosity 
and the one based on the resulting strain rates. With MARC this turned out not 


to be necessary except for n < 1. 
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To reduce computation time the effective viscosity of every integration point of 
every element was written to a file and during subsequent calculations the resulting 
effective viscosity of the most similar calculation was used as a starting value, which 
led to substantial reduction of CPU-time. 


6.3. Mesh generation 


Since a large number of models with various different € and 6 values was to be calcu- 
lated, it was impractical to generate a new FE mesh for each and every calculation. 
Instead a small number of rectangular meshes were generated, which were then — 
through the use of transformation formulae — transformed into the desired shape. 


The transformation formulae were 


zg =ae~*/h-ek/2 (6, 4 & — Ee? — foe**) sin kx + 2, (6.1) 
z =a(te** — &e**) coska + (1+ ak( + &)(e7/*-e-"/*), (6.2) 

where ohh 
= ah and 9 := -14+&, (6.3) 


(x ,z) denote the transformed values and (z, z) the original ones. Note that the 
bed was defined by 
2 = acoskr (6.4) 


and not z) = asinkz. In some cases the mesh was scaled linearly in both directions 
prior to the use of (6.1) and (6.2), keeping the wavelength fixed, to obtain different 
wavelength-to-thickness ratios. In this way it was also possible to apply different 
scalings to get different meshes with the same ¢ and 6 values, making it relatively 
easy to estimate the numerical error associated with a change in mesh geometry. 
The transformation formula was found by starting with the formulation: 


z = f(z)acoskz + g(z) (6.5) 


where f(z) and g(z) are some unknown functions satisfying: f(0) = 1, f(h) = 0, 
g(0) = 0 and g(h) = 1. It is desirable that the transformation formula gives ap- 
proximately a conformal mapping, but it should not be an exact conformal mapping 
because there will inevitably be a “crowding” of the transformed 2'-lines, giving 
bad aspect ratios for the generated elements (Menikoff and Zemach, 1980). The 
Cauchy-Riemann equations are 
zt, = z and L,= —z, (6.6) 
Putting (6.5) into the left-hand side of (6.6) and integrating gives 
t =akF(z)sinkz + K(z), (6.7) 


where F is the integration of f and K is an unknown function. Using the right-hand 
side of (6.6) gives a differential equation for f(z), which can easily be solved 


f(z) = aye"? — age™, (6.8) 
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where a, and ay are constants, and g. = A, which gives g(z) = a3z and A (rc) = a3z. 
It is required that for z = 0, z is transformed to zo, and to z =hforz =h. In 
order to avoid the negative aspects of crowding the factor e~*/4~@*/? was introduced, 
which makes the transformation non-conformal. It then turned out to be necessary 
to ensure that the transformation is one-to-one’, which was guaranteed by requiring 
z, > 0 and x, > 0. Eq. (6.1) and Eq. (6.2) can be seen as a compromise between 
the positive aspect of a conformal mapping giving small skewness values for the 
generated elements, and the negative aspect of the crowding, which gives large 
aspect ratios. This transformation gave an excellent FE mesh and could be used for 


all values of € and 6 of interest. 


6.4 Boundary conditions 


The top of the model is a free surface; there o,, and o,, are equal to zero. The 
atmospheric pressure was ignored. 


The perfect sliding condition was applied to the bottom nodes. It was implemented 
by suppressing all movements in the direction normal to the bed and allowing free 
movements in the tangential direction. 


At the sides a periodic velocity boundary conditions was specified, 7.e. the left side 
was forced to have the same velocity profile as the right side. The model is infinite 
and periodic. No prior knowledge of the velocity or of the stress tensor along the 
edges was needed. 


6.5 Post-processing 


There were basically two types of post-processings performed. One was the direct 
visual inspection of the results of an individual calculation with the post processing 
program Mentat. This was time consuming and could only be done for a limited 
number of runs. The other one consisted of an automatic calculation of all quantities 
of interest. Various subroutines from the IMSL collection of numerical algorithms 


were applied in doing so(IMSL, 1989). 


A one dimensional periodic B-spline interpolation through the velocity values at the 
bed was made to obtain equally-spaced velocity values along the x axis, which were 
then Fourier transformed and used for calculating the mean of v,, te. the sliding 
velocity. 


The velocity field was interpolated by the use of a bilinear interpolation function 


that was consistent with the FE mesh and the FE interpolation functions used. 


The searching procedure for Une and Us", involved a non-linear constrained min- 


imization/maximization problem because the = extension of the region of interest 
. . see . rmax , 
is dependent on the value of x. As starting values for the positions of U7" and 


1A function T: V + W satisfying T(z) = T(y) and implying r = y, is said to be one-to-one 


on VY. 
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Us the analytical results discussed in Sec. 4.3.2 were used. The positions of Une 
and aa within the FE mesh were also found using the nodal velocity values and 
compared with those found by interpolation. The difference was never large, which 
is the expected results, or comparable to the size of the elements. By interpolation 
of the nodal values somewhat more precise results were obtained than otherwise, as 

no and aD could then be found somewhere within the elements. This allowed 
a better comparison between theoretical and numerical results. 


6.6 Numerical results 


6.6.1 Sliding velocity 


In Sec. 3.2 it was shown that in the limit as ¢ > 0 and for 6 < 1 the sliding function 
s, defined by Eq. (3.24) will be a constant co that is independent of ¢ but depends 
on n, 1.€. 

lim s(e, 6, n) ‘xl Co(d, 7), (6.9) 
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and it was argued that co = ¢'~” with ¢ = 2. Otherwise no estimates of the sliding 
function s could be made, and it was realized that numerical work could give some 
more information about s. 


The sliding function was calculated for n = 1 to n = 5, for r = 0.001 tor = 10 
and for ¢ in the range from 0.025 to 1.0. Some of the results can be seen in Fig. 6.1 
and in Fig. 6.2, which show Ins as functions of Ine for n = [1,2,3,4,5] and for 
6 = 0.05/(27) « 1 and 6 = 0.1/(27) respectively. The figures display a number of 
interesting features: 


e For small ¢ values all curves have zero slope. Eq. (6.9) is therefore correct, and 
the sliding velocity is, for € small enough, proportional to e~('+™ as Eq. (3.23) 
indicates. This can be seen more clearly in Fig. 6.5-a that shows the function 
U,, or the nondimensional sliding velocity (defined through Eq. (3.6)), for 
€ < 0.125.? The slope of the f(e) curves is given in Table 6.2, and depicted 
as a function of n in Fig. 6.5-b. The calculated slope is —1.017 — 0.986 n, or 
within 2 % from the theoretical slope of —(1 + n). 


e The range of the validity of Eq. (3.23) depends on n. As n increases this range 
of € values gets smaller, contrary to statements made by Fowler (1981, p. 675) 
who argued that e”*! had to be small. For n = 5, ¢ must for example be smaller 
than about e~? = 0.135 (r < 0.0215) for Eq. (3.23) to be a approximation valid 
within 20%. Although a value less than 0.02 for r is not unreasonably small, 
this fact considerably reduces the usefulness of Eq. (3.23). For n = 1 the 
variation of s is slow; not until ¢ = 1.0 is s 20% different from s(0, 0.05/27, 1). 


For € and 6 fixed Ins changes by a constant amount for every unit change in 
n, te. In(s(e,d,n + 1)) — In(s(e,6,n)) = A(e, 4). It is an important fact that 


ee St te a ee 
?Uy and s are connected through s = me"+! Uy. 
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Fig. 6.1: The logarithm of the sliding function s as a function of Ine and n for > = 0.05. 
Every symbol represents the result of one calculation. The constant slope of the curves for 
€ —+ 0 shows that u, o« e~("+) in that limit. The range of « values for which this asymptotic 
expression is correct depends on n, and becomes smaller with increasing n. Note that for every 
particular value of € there is a constant shift in the value of Ins with each unit increase in n. 
This means that Ins is a linear function of n. It is therefore not necessary to calculate s for 
every possible value of n; all that is needed is a simple interpolation or extrapolation of the 
calculated values. 


A does not depend on n. Hence, one does not have to calculate s for every 
possible value of n. A simple interpolation or extrapolation based on several 
different n values will suffice. In Fig. 6.3 the logarithm of the sliding function 
s is shown as function of n for a few roughness values r and a fixed corrugation 
value of ¢ = 0.005. It is seen that a straight line always results, but with a 
slope dependent on ¢. For ¢ = 0.005 (the smallest value of ¢ in the figure) 


In 5(0.005, 0.05/27, n) = 0.5138 — 0.5370n (6.10) 


with a standard deviation of 0.010. Eq. (3.29) predicted this power-law be- 
havior but with a slope of —In2 ~ —0.69. Notice that relation (6.10) gives 
Ins = —0.0232 4 0 for n = 1, although s was defined in such a way that 
In s(0,0,1) = 0. This discrepancy is not entirely a result of numerical errors 
since it was obtained with finite values of ¢ and ¢. 


For 6 = 0.05/27, s was approximated as 
6 . 
s(e,d,n) = do ca (6, n)e” (6.11) 
1=0 


for ¢ < 1/2. The resulting curves are seen in Fig. 6.4, which shows s as a 
function of <2, and the values of the Taylor coefficient are given in Table 6.1. 
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Fig. 6.2: Ins as a function of Ine and n for ¢ = 0.1. The only difference between this figure 
and the previous one is the value used for 6. No difference can be seen, showing that the 
results do not depend on the exact value of 6. 


In(s) 





Fig. 6.3: Ins as a function of n for several roughness values r and ¢ = 0.05. For every 
roughness value r a linear relationship between In s and n is found, with a slope that depends 


on r. s is defined in such a way that the limiting value of Ins for 7 3 0 and n = 1 is equal 
to unity. 


If s is needed for some other n, s should be calculated according to Eq. (6.11) 
and interpolated or extrapolated using the fact that In s(,,n) = é +¢.n where 
¢, and ¢p are obtained through a least squares approximation. 


Attempts to extrapolate to ¢ = 0 did not convincingly result in improvements, 
partly because practically the same slopes were obtained for different 6 val- 
ues as long as 6 < 1 but more importantly because s was calculated for a 
rather limited number (+ 10) of wavelength-to-thickness ratios. Comparison 
of Fig. 6.1 and Fig. 6.2 shows that a factor of two change in 6 does not change 
s significantly, 2.e. € is much more important for s than 6 (How s depends on 6 
will be discussed in Sec. 6.6.2.). By calculating the slope of the s(e) curve for 
different values of 6, it could be shown that the small but finite value of 6 is 
responsible for the 0.64 % deviation of co(0.05/27, 1) from 1 seen in Table 6.1. 
Although this deviation is small, a standard statistical test showed it to be 
significant. 

The values of co for n = 3 from Table 6.1 is in agreement with the esti- 
mate (3.26) from Meyssonnier based on theoretical arguments. The value of 
c2(0,3) = 1.163 + 2% from Table 6.1 is, however, most probably not in ac- 
cordance with his estimate for c2(0,3) = 2.4. This deviation could be caused 
by the limited number of calculations, numerical errors, or both. Despite this 
difference there seems to be no reason to doubt the general correctness of 
Meyssonnier’s findings. 


The numerical results are also in agreement with Lliboutry’s (1993) theoretical 
upper estimate for s(0, 3) 


s(0,3) < 0.33839 + 3.688? + 0.17 €7 (6.12) 


and his tabulated values for s (Lliboutry, 1993, p. 62). His upper estimate is, 
however, up to several times greater than the calculated values.? 


All curves in Fig. 6.2 intersect at Ine ~ —0.3. This fact does not need any 
special attention. If one wants to compare s for different n as a function of €, 
s would have to be made comparable by, for example, normalising it with cp. 
Then all one could say is that s(e,6,n +1) > s(e,6,n). 


For Ine > 1, U, x €?, where p depends on n. This can be seen in Fig. 6.5- 
c, and in Fig. 6.5-d where 0InU;/Olne for Ine > 1 as a function of n is 
depicted. The best straight line approximation is given by —1.11 — 0.228n 
with a standard deviation of 0.010, i.e. Uy « ¢~':!!~°-228"_ Hence, for extreme 
roughness the sliding velocity is not dependent in the same way on € as it is 
for slight roughness. 


s(e,6,n) as a function of € and ¢ for n = 3 is shown in Fig. 6.6. Diamond symbols 
represent the (¢,<¢)-pairs for which s was calculated. The lines were calculated by 
interpolating this randomly scattered data. It can be seen that s depends more 
strongly on ¢ than ¢, and that 0s/0¢ is negative and increases in magnitude as ¢ 


increases. 


3Lliboutry (1993) uses the symbol V’ for the sliding function. 
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Fig. 6.4: s as a function of <? for 6 = 0.05/20 and « < m/2 (7 < 0.25). The symbols 
represent calculated values and the lines are least squares approximations using s(e,6,n) = 
37° 4 Coi(6,n) €7*. Table 6.1 gives the values for c9;. 


pel 4 5 
0.9936 |] 0.5661 | 0.3294 | 0.1943 
-0.1486 | 0.9328 | 1.163 | 1.055 
0.1137 | -0.4997 | -0.3795 | 0.06465 
0.08996 | -0.3361 | 0.3066 | 0.4760 


-0.01813 | 0.5192 | -0.4693 | -0.7965 
-0.001253 | -0.2334 | 0.2459 0.3642 

0.0006450 | 0.03723 | -0.04077 | -0.05325 
3.5 x 107° | 0.0028 | 0.0027 0.0067 








Table 6.1: Taylor coefficients of the sliding function s(e,6,n) = oo cai(6,n)e” for 6 = 
0.05/27. These values can be used for 5 = 0 with less than 2% error. 


n Qn bn standard deviation 
1 -1.159 -2.002 0.00071 
2 -1.671 -2.987 0.0017 
3. -2.174 -3.976 0.0032 
4 -2.653 -4.962 0.0047 
5 -3.110 -5.943 0.0060 


Table 6.2: Linear regression coefficients for InU,(¢,0.05/27,n) = an + ban for ¢ < 0.125. 
Theory predicts b, = —(n + 1) for e 3 0. 





Fig. 6.5-a: InU, as a function of ¢ and n for 
¢ = 0.05 and € < 0.125. 

The straight lines show the best linear 
approximations through calculated values 
given by the symbols. 
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Fig. 6.5-c: InU, as a function of e and n for 
¢ = 0.05 and € > 2.7. 

The straight lines show the best linear ap- 
proximations though calculated values given 
by the symbols. 
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Fig. 6.5-b: U, as a function of n for ¢ = 0.05 
and € < 0.125. 

Each point represents the slope of a curve 
in Fig. 6.5-a for the corresponding n. The 
best straight-line approximation is given by 
—1.017 — 0.986n with a standard deviation 
of 0.0025, in an agreement with the theo- 
retical prediction of —(1+ 7) for the slopes 
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Fig. 6.5-d: U, as a function of n for ¢ > 2.7 
Each point represents the slope of a curve 
in Fig. 6.5-c for the corresponding n. The 
best straight line approximation is given by 
—1.11—0.228n with a standard deviation of 
0.010. No theoretical prediction known. 


The sensitivity of the sliding velocity to changes in glacier thickness is therefore 
larger for high corrugation values than for small ones. 


6.6.2 The dependency of wu, on 6 


The contour lines in Fig. 6.6 are more or less straight lines (deviation from straight 
lines could be a result of the interpolation procedure) and this fact indicates that s 
is possibly a linear function of 6. By plotting s as a function of ¢ for different values 
of n and ¢, as done in Fig. 6.7, it is found that s is indeed a linear function of ¢ and 
that the slope increases with increasing ¢ and, more importantly, with increasing n. 
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Fig. 6.6: Dependency of s on € and ¢ for n = 3. Contour interval is 0.05. Diamond symbols 
represent the (€,¢) pairs for which s was calculated. It can clearly be seen that s depends more 
strongly on « than ¢, and that Os/0c is negative and increases in magnitude as € increases. 
As explained in the text, a closer inspection of the numerical results showed that the sliding 
function s depends linearly on 6. The contour lines should therefore be straight lines, and the 
deviation from that form results from errors associated with the interpolation. 


Hence, the sensitivity of the sliding velocity to changes in glacier thickness increases 
with increasing n. Except for n = 1, s is a decreasing function of 6 and wu, therefore 
decreases somewhat faster than 7;’. 


For r = 0.02, s changes less than 6% for 0 < ¢ < 1 even for n = 5. Fors & 1 the 
variation of the sliding velocity is expected to set up some surface variation that 
in turn will change the stress pattern at the bottom resulting in a different sliding 
velocity. ¢ < 1 can therefore be seen as a part of the definition of the problem 
and not a simplifying assumption. If ¢ < 1 is not fulfilled the notation of sliding 
velocity is not clear. For reasonable ¢ values (¢ < 1/4) ¢ does not influence the 
sliding velocity significantly unless n > 3 and r > 0.1. As an example, for n = 5 and 
r = 0.2, s changes ~ 20% as ¢ goes from zero to 0.25. This is unexpectedly large, 
and it is important to get some rough idea of what is causing this ¢-dependency. 


The following simple explanation turns out to be helpful. For ¢ < 1 the contribu- 
tions of the crest (0 < X < 7) and the trough (m < X < 27) of the sine wave to the 
average shear stress are equal. As € increases this is no longer true as higher har- 
monics that disturb this symmetry start to be important. The stresses become miore 
and more concentrated about the crest meaning that the force balancing gravity is 
increasingly supported by it. The glacier thickness, as well as the driving force, are 
therefore effectively reduced. Let us assume that the driving stress is given by 


T = pglh — p(r)a) sina (6.13) 
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Fig. 6.7: s as a function of ¢ for r = 0.02 (solid lines) and r = 0.2 (long dashes). The numerical 
results indicate that the sliding function s depends linearly on the \/h ratio (asperity). The 
slope of the s(¢) curve depends on n and e, and it can be shown that this dependency is 
approximately linear, and that s(e,6,n) (1 — e6n)s(e,0,n), for 6d < 1. 


with 0 < p(n) < 1, where p(x) is otherwise an unknown function, which depends 
on n because the shift in stress pattern and the degree of stress concentration is 
expected to depend on n. If the only influence of a change in a is given by this 
reduction in driving stress then it follows that 


s(e,¢,n)/s(e,0,n) = (7/%)” 
(1 — edp(n))” (6.14) 


@S! (1 — ednp(n)) 


The slope of the s(5) graph should therefore be proportional to ¢ and nr(n). The 
ratio s/(1 — €d)" as a function of ¢ is seen in Fig. 6.8. Comparison with Fig. 6.7 
shows that the slopes are greatly reduced, meaning that the most part of the 6 
dependence can be explained in this way with p(n) = 1. The systematic change of 
s/(1 — €5)" as a function of the asperity ¢ is (for a reasonable range of ¢) so small 
that it is of no practical importance. The curves, however, seem to be somewhat 
over-corrected for small n values, and under-corrected for large ones. This would 
mean that p(n) increases with n or that the stress concentration becomes larger for 
increasing n. An inspection of the numerical results showed this to be true. 


Of interest is the sliding velocity as a fraction of the total surface velocity. If one 
writes the surface velocity u, as a sum of the deformation velocity ua and the sliding 
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Fig. 6.8: s/(1 — £6)” as a function of ¢ for r = 0.02 (solid lines) and r = 0.2 (long dashes). 
By comparing this figure with the previous one it can be seen how much of the dependency of 
s on 6 can be accounted for by the assumption s(e, 6,n) = (1 — €6)"s(e,0,n). 





velocity* 
ai? s 
Us = 2A Tp (5 + =) (6.15) 
one finds 
up ett] a1 6.16 
Ay := — = {1+ ———_ ] . ; 
5 uy ( - (n+ 55 i) 
The ratio that determines the fractional sliding velocity (A,) is therefore: 
entl 
a 6.17 
(n+ 1)6s’ cen) 


and A, can be quite large. For r = 0.1, 5 = 0.05 and n = 3 one finds for example, 
using the values from Table 6.1, that A, = 0.58. Using r = 0.05 gives A, = 0.99. 
Sliding over a hard bed without bed separation can be significant. It is important 
to realize that Eq. (6.16) gives the maximum fractional sliding possible. In nature 


ice will always experience some friction from the bed reducing the sliding velocity 
(possibly down to zero). 


Eq. (6.16) shows that 0A,/06 depends on €. Observation of changes in fractional 


sliding velocity as a glacier changes its thickness could therefore give information on 
the roughness of the glacier bed. 





‘This is only approximately correct. For n # 1 the ice close to the bedrock becomes somewhat 
softer and the internal deformation therefore somewhat. larger than otherwise expected. 
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Fig. 6.9: The velocity increase of U7"3* with respect to the velocity at bed at X = 1/2, for 
¢ = 0.05. With increasing «, U™3* moves towards Use and disappears at € = Ecriticay. For 
n= 1 and 6 < 1 this limiting value is approximately equal to 0.138. The vertical dashed- 
dotted line marks this « value. For n = 1 there should therefore theoretically be no symbols 
to the right of the line. The numerical calculations confirm this. The value of ¢ above which 
no local velocity maximum was found depends on n. This is reflected in the figure by the fact 
that symbols are found at progressively larger « values with increasing n, showing that the 
magnitude of the velocity maximum is both enhanced by the non-linearity of the flow law and 
that it exists up to larger roughnesses. 


6.6.3 Extrusion flow and non-linear material behavior 


As explained in Sec. 4.3 extrusion flow is expected to occur above the point where 
Um for n = 1 arises, as long as € < 0.1378839, and below the point where ae 
arises as long as € < 1/2. Extrusion flow is caused by the bedrock undulations 
and can therefore be seen as a geometrical effect. The flow law can, however, be 


expected to have an influence on the magnitude and the onset of extrusion flow. 


In Fig. 6.9 the percentage of increase of the velocity U7"5* over the velocity at Zo(1/2) 
as a function of ¢ is shown for n = 1 ton =5. The vertical dashes-dotted line is at 
€ = 0.1378839. To the right of that line no plus symbols are shown, indicating that 
in agreement with the theoretical prediction for n = 1 no U7"3*-points were found 
for larger ¢ values. The long dashes show the theoretical asymptotic slope for « — 0. 


Again a good agreement between numerical and analytical results is found. 


For n > 1 extrusion flow exists. As n increases the range of ¢ values for which 
extrusion flow is found becomes progressively larger. The velocities also become 
progressively larger. Extrusion flow is therefore enhanced by the non-linearity of 
Glen's flow law. For n = 3 it is for example found for € up to © 0.2 or r = 0.03, and 
the velocity is © 10% larger than at the bed. The position of U7/3* for n = 1 follows 
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Fig. 6.10: Relative decrease of Oss with respect to the velocity at the bed at X = 3n/2 
for ¢ = 0.05. For ¢ — 0 there is no velocity variation and the velocity of the local velocity 
minimum is the same as the velocity at bed. Thus for « — 0 all the curves approach the value 
0. For € > 0 the minimum velocity is always less than the velocity at the bed (extrusion flow) 
so that the percentage increase is always negative. With increasing e, Us becomes, as a 
percentage of the bed velocity, larger negative, and at some particular value of € a minimum is 
reached. The position of this minimum moves to larger ¢ values with increasing n. The velocity 
minimum is also stronger for larger n values (larger negative numbers). The non-linearity of the 
flow law therefore has the effect of enhancing the vertical velocity variation in relative terms, 
and with it the extrusion flow behavior. On the basis of Morland’s solution it was expected 
that for a linear medium no velocity minimum would be found for « > 1/2. The results of the 
numerical calculations do not, as can be seen in the figure, support this conclusion. Up's for 
n = Lis found up to e = 1.18. The reason for this discrepancy lies presumably in the fact that 
1/2 & 1, and that therefore the results of the second order perturbation theory, were ¢ < 1 
is assumed, cannot be trusted. 


closely the theoretical prediction. For n > 1 the 7/2 Moves progressively closer to 
bed. This is in agreement with other numerical calculations (Raymond, 1978). 


Fig. 6.10 shows the minimum of the absolute horizontal velocity at X = 31/2 as a 
percentage of the velocity at the bed 2.e. 


100 (Meret ba 7 (6.18) 


for ¢ = 0.05 as a function of ¢. The theoretical prediction for n = 1 was that Up 
should exist for ¢ < 1/2. Symbols below the horizontal zero line indicate that a 
minimum was found. For n = 1 a minimum was found for ¢ up to 1.17 and not. only 
up to ¢ = 1/2. The theoretical value is based on a perturbation analysis which is 
only valid fore < 1. Since Ut exists for ¢ > 1 it is clear that the perturbation 
approach could never have given the correct answer. The asymptotic change of the 
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Fig. 6.11: Vertical position of U. Use (s = 0.05). The numbers on the y-axis correspond to a 
vertical distance of 1/k. For « + 0, Up i is situated at Z = 1, or z = 1/k. If the minimum 
velocity at kx = 37/2 is found at z = —a, UR is at the bed. This corresponds to kz = —ka 
or Z = —e. The line made of points and long dashes gives the position of the bed as a function 
of e. If the symbols lie on the line there is no local velocity minimum above the trough of the 
sine curve. Note that with increasing n, Uj7'7, moves away from the bed-tine. 


velocity decrease as ¢« — 0, shown as a solid line, is however reproduced. There 
is therefore a good agreement between theory and numerics at € values where an 
agreement can be expected. Notice that the two figures are practically identical 
indicating that the exact value of ¢ is unimportant as long as it is small. This is 
also correct for others figures shown below. 


For all calculated values of n the velocity decrease (shown in the figures as a nega- 
tive velocity increase) becomes larger as ¢ increases from zero, reaches a maximum 
and decreases again. There is always some ¢ value above which no extrusion flow 
is found. As n increases the extrusion flow is strongly enhanced. It again be- 
comes progressively larger in magnitude and exists for larger € values as n increases. 
Comparison of Figs. 6.10 with Fig. 6.9 shows that extrusion flow behavior is more 
dominant above a trough than above a riegel. Above a trough a 30 to 40 % decrease 
in horizontal velocity with height over a distance of ~ 1/k is possible and could for 
example cause a considerable inversion of a bore-hole inclination. 


U. tay arises for ¢ = 0 at Z = 1 and moves towards the bed with increasing ¢. How 


this happens for n = 1 to nm = 5 and ¢ = 0.05 can be seen in Fig. 6.11. The 


dashed-dotted line gives Z = —e, which is the position of the bed. As long as the 
symbols remain above the dashed-dotted line (3 exists. It is interesting that U ery 
is situated slightly larger above the bed for larger n values which is opposite to the 
behavior of U7/5". It is therefore not a general rule that characteristics of the ice 


deformation iad to approach the bed with increasing n. Another interesting feature 
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Fig. 6.12: Velocity at (X, Z) = (37/2, —e) asa fraction of the sliding velocity (¢ = 0.05). Note 
that the velocity at the base of the trough at first increases with ¢, and that the non-linearity 
of the flow law makes this increase larger and displaces it to larger € values . 


of the figure is that Say remains more or less at a constant height of approximately 
1 (or at 4/27 in dimensional units) above the bed for all values of e. 


The velocity within an overdeepening as a fraction of the sliding velocity is often of 
interest. It is for example sometimes important to know at what roughness values 
the ice within the overdeepening effectively remains there without taking part in the 
overall glacier motion. What exactly is meant by saying that the ice does not move 
depends of course on what part of the overdeepening — at the bed or only “close” 
to the bed — one is referring to, but it turns out that this is, at least for the case 
of a perfectly lubricated bed, relatively unimportant, as will be shown below. 


The velocity at (X, Z) = (32/2, —e) as a fraction of the sliding velocity (cf. Fig. 6.12) 
becomes, as € increases from zero, at the beginning somewhat larger. This unex- 
pected result is a manifestation of the extrusion flow and is more important for 
non-linear than for linear flow behavior. As € increases further the ratio v(X = 
3n/2,Z = —e€)/us decreases fast; it is less than 0.1 at ¢ = 1.6 and negligibly small 
for € > 2. 


Another measure of the magnitude of ice movements within a trough would be 
the velocity at the base of the trough ((X,Z) = (37/2,-e)), or the minimum 
velocity at X = 32/2 for some Z, normalized by the velocity at the top of the riegel 
((X,Z) = (m/2,¢)), which is depicted in Fig. 6.13 and Fig. 6.14 respectively. In 
all cases the result is effectively the same. For ¢ > 2 there is practically no ice 
movement within the trough. 


A value of ¢ = 2 corresponds to r = 1/m & 0.32 which can hardly be considered to 
be a very high roughness value. For an overdeepening or a trough having this value 
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u(X=3n/2,Z=-e)/u(X=n/2,Z=.) 





Fig. 6.13: u(X = 30/2, Z = —e) normalized by u(X = 7/2, Z =e) for ¢ = 0.05. The curve 
gives the velocity at the base of the trough as a fraction of the velocity at the peak of the sine 
curve. 


of roughness or larger the ice within it practically does not move. Perfect sliding 
is an idealization, which will cause the maximum local basal velocity possible. Any 
degree of friction will presumably pose a further hindrance to ice movements. 


The minimum velocity at X = 37/2 normalized by wu, is shown in Fig. 6.15. For 
€ large enough to exclude extrusion flow the minimum is found at the bed (Z = 
~—e), and Fig. 6.15 shows a similar behavior Figs. 6.12 and 6.14. But there is also 
something quite unexpected to be seen; for € > 1.8 u(X = 3n/2,Z = —e)/us is 
negative, i.e. the medium in the lowest part of the trough flows in the opposite 
direction to the main flow. This is a clear indication of a flow separation. 


An example of a flow separation is given in Fig. 6.17. It is an enlarged part of 
Fig. 6.16, which shows the flow above and within an overdeepening for A = 50m, a = 
20m, h = 200m, n = 1 and pgsina = 8.99577 x 10-* bar/m. Fig. 6.16 again only 
shows a part of the hole configuration, which had the dimensions 200 x 200m. The 
velocities have the dimension m/a. Although the velocities within the overdeepening 
are small compared to the velocity at the riegel, they are a significant fraction of 
the sliding velocity, which for this particular case is 0.50 m/a. It can clearly be seen 
how the main flow induces a secondary flow circulation in the clockwise direction. a 
separation line is formed (shown as long dashes), separating the main flow from the 
induced flow. The ice below the separation line will theoretically circulate there for 
ever, never leaving the trough. 


Flow separation is not, as sometimes thought, limited to high Reynolds numbers as 
Figs. 6.17 and 6.16 clearly show. There is also a number of analytical and experi- 
mental studies of creeping flow that show that flow separation is not limited to the 
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min, (u(X=32/2,Z))/u(X=1/2,2=e) 





Fig. 6.14: The minimum of the horizontal velocity v, at X = 37/2 for some z normalized by 
the velocity at the peak of the sine curve, i.e. minz(u(X = 37/2, Z))/u(X = 71/2,Z =e). 
The minimum is not necessarily found at the base of the trough. The negative numbers 
indicate that the minimum velocities are negative, i.e. the flow direction is opposite to the 
flow direction at (X,Z) = (7/2,¢). This is an example of the recirculation of the ice within 
the trough which takes place at high < values. 


particular case of a perfectly lubricated sinusoidal bed (Hasimoto and Sano, 1980; 
Sherman, 1990, p. 258 - 265). Corner flow, as an example, driven by circumferential 
motion with no-slip boundary conditions, is known to form so called Moffatt corner 
eddies® if the corner is sharp enough (Moffatt, 1964). The corner angle must be less 
than 73°. A review of the literature has, however, not revealed any other examples 
of flow separation for perfectly-sliding-type boundary conditions. Separation is also 
known to occur in creeping flow past sharp-cornered obstacles such as a vertical 
fence of in a plane flow past a step (Taneda, 1956). 


Tsangaris and Potamitis (1985) did numerical calculations on laminar small Reynolds- 
number flow over a sinusoidal wall using no-slip boundary conditions. They found 

flow separation at ¢ ~ 0.6 for Re = 1, where the Reynolds number Re was defined 

by Re = UX/(2n), where U was the horizontal velocity at large distance above the 

wall. 


The magnitude of the reversed flow can be several times as large as the sliding 
velocity (cf. Fig. 6.15), but enlarging n decreases it. To test the influence of the 
FE grid size on the onset of flow separation a local mesh refinement was done that 
made the mesh unsymmetrical with respect to X = 37/2. It was found that this did 
locally increase the reverse flow. The possibility that extrusion flow occurs somewhat 
earlier than Fig. 6.15 indicates and that its magnitude is somewhat greater cannot 


°Named after H. K. Moffatt, who theoretically proved their existence 
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min, (u(3x/2,2))/u, 
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Fig. 6.15: The minimum of the horizontal velocity above the trough of the sine wave as a 
function of € (€ := ak), ie. minz(u(3m/2,Z)). The velocities were normalized by the average 
sliding velocity u,. For ¢ — 0 there will be no velocity variations and the minimum of v, is 
equal to uy. Therefore all the velocity curves start at 1.0 for € = 0. With increasing roughness 
the minimum of the horizontal velocity as a fraction of the sliding velocity becomes smaller. For 
€ = 1.75 the minimum of flow velocities within the trough have become a negligible fraction 
of the sliding velocity. As € increases further negative values are found. This shows that the 
minimum velocities have become negative. Since the roughness is then quite large the sliding 
velocities are small and the ratio of the minimum velocities (maximum negative velocities) to 
uy, can thus be targe. 


be excluded. This is especially true for n > 1. 


6.6.4 Frequency doubling 


That larger frequencies become increasingly important can clearly be seen in Fig. 6.18 
and 6.19 that show the first four harmonics of the basal velocity for n = 1 and n = 3. 
The figures are double logarithmic and the Fourier coefficients have been normal- 
ized by uy. The solid line gives the increase of the second harmonics according to 


Kg. (4.8). 


For n = 1 (Fig. 6.18) and € small, the first two harmonics are about equally impor- 
tant. But as € increases and the stresses concentrate about the peaks of the sine 
curve the amplitude of the other two harmonics shown become comparable with 
the first two ones. Within the framework of perturbation analysis each harmonic is 
associated with one term of the power series. The rise of the i-th harmonic with € 
is therefore proportional to the i-th power of €. This can be seen in Fig. 6.18 where 
the slope of the curves is equal to the number of the harmonic. 
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Fig. 6.16: Flow above and within an over-deepening. Parameters used: 4 = 50m, a = 20m, 
h = 200m, n = 1 and pgsina = 8.99577 x 10~? bar/m. The velocities have the dimension 


m/a. Only a part of the FE model is shown. 


For n = 3 the higher harmonics turn out to be relatively more important with 
respect to the first harmonic than for n = 1. Frequency doubling seems to play a 
more important role for non-linear flow. This can be understood as a result of the 
flow becoming more concentrated about the peaks for n > 1. 
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Fig. 6.17: Flow separation. Same parameters used as in the previous figure. Only part of the 
FE mesh is shown. The direction of the main flow is from left to right. The vectors indicate 
the direction of the flow at each FE node. The ice within the trough rotates slowly in the 
clockwise direction. The ice which takes part in this circular motion will theoretically never 
leave the trough, and is thus separated from the main flow. This calculation was done for a 
perfectly lubricated bed. No systematic study was made of the effect of the boundary condition 
along the bed-line on the flow behavior. This particular calculation was, however, repeated for 
no-slip boundary conditions, and again a flow separation of this type was found. 
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Fig. 6.18: Harmonics of the basal velocity for m = 1 (¢ = 0.05). The Fourier coefficients 
were normalized by the sliding velocity u,. The solid line gives the theoretically-derived curve 
for the second harmonics. The curve is double logarithmic. Note how important the second 
harmonics are. With increasing € all harmonics become equally important. 
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Fig. 6.19: Harmonics for n = 3 (¢ = 0.05) normalized by u,. The solid line gives the 
theoretically-derived curve for the second harmonics for n = 1. The curve is double logarithmic. 
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CHAPTER 7 


Conclusions 


With the use of both analytical and numerical methods the flow characteristics as 
well as the sliding velocity of a highly viscous medium flowing under the influence 
of gravity over a perfectly lubricated sinusoidal bed have been analysed. In this 
chapter a summary of results, and recommendations for further research will be 
given. Examples of field observations of extrusion flow that can be understood in 
the light of these new findings and predictions about flow patterns of glacier in the 
Alps that could not have been done before this investigation was made, will be given. 


7.1. Flow characteristics 


The horizontal velocity field of a highly viscous medium flowing over a frictionless 
sinusoidal bed can have a point of local velocity maximum 7/2 above the peak 
of the sine curve, and a point of local velocity minimum above the trough of the 
sine curve. Whether these stationary points will exist or not depends mainly on the 
roughness of the bed (r := a/A or € := ak), and the non-linearity of the medium, 
but also to some extent on the wavelength to thickness ratio. 


7.1.1 The local velocity maximum, a2. 


The velocity maximum is found for € < Ecritica(1™/2, 5). Ecritical is a function of n and 
is listed in Table 7.1 for n = 1 ton = 5 ford < 1. The value of €citicay increases 
with increasing n ppowine that the non-linearity of the flow law makes the range of 


€ values for which U7"5* exists larger. 


n/2, 8 Situated directly above the peak of the sine curve. For € = 0 its vertical 
position is given by Z = 1. With increasing e, U? r/2 moves upwards with respect to 
the mean bed-line. For n = 1 and € = Ecritical, Uz/2° is at Z ~ 1.98. Increasing the 
value of n moves Un towards the bed. 


There is a point of horizontal velocity field situated directly above U7/5" were v; 
has a local maximum in the horizontal direction but a local minimum in the vertical 
direction. This saddle point Us" moves downwards from Z = 00 towards the 
maximum point Un with ificroueing e. For € = Ecritical the vertical position of the 
saddle point and that of the maximum point become coincident. For € > Ecritical, 


mr and US Loria do not exist. 


a9 A 


nr Ecritical(™/2, 0) 
1 0.13788 

2 0.19 

3 0.22 

4 0.24 

5 0.25 


Table 7.1: €criticat(77/2,0) as a function of n for d < 1. The value of Eciticar for n = 1 is 
based on the Morland solution. The values for n > 1 are based on numerical calculations. The 
velocity maximum above the peak of the sine curve (ka = 7/2) exists for € < Ecritical. The 
errors in the determination of Ecriticat are estimated to be less than 7%. 


Between Usage and U7" the horizontal velocity increases with depth. This is an 
example of extrusion flow. Extrusion flow would cause an inversion of the inclina- 
tion of bore holes. It is interesting that this type of inclination inversion has been 
observed in nature. Inclinometry measurements in a bore hole slightly down-glacier 
from the crest of a riegel on Storglaciaren, Sweden, showed an increase in horizontal 
velocity with depth (Hooke et al., 1987). This was only observed during summer 
when water pressures were high. There are indications that the ice became at least 
partly decoupled from the bed during this period. The assumption of frictionless 
flow over the bedrock may thus be valid. A map of the bed topography obtained 
by radio-echo measurements (Bjérnsson, 1981) shows that there are two regions of 
overdeepening up-glacier and down-glacier from the riegel. The bed can thus be 
approximated by a sine curve. The wavelength is about 1000 m and the amplitude 
approximately 20 m which gives ¢ + 0.12. This € value is less than €criticai(1/2, 0) 
from Table 7.1 for all values of n. The assumption of 6 < 1 is however not valid and 
one can therefore not expect a perfect numerical agreement. The measured vertical 
position of the region of increase in horizontal velocity with depth is, for example, 
closer to the bed than would be expected on the basis of the numerical calculations 
done for 6 < 1. Still, it can be concluded that the vertical extension and compres- 
sion of the ice caused by the sinusoidal form of the bed is a likely explanation of the 
observed flow pattern. 


7.1.2 The local velocity minimum, a) 


The minimum point aero is situated above the trough of the sine curve (kx = 37/2), 
and is found for € < €critical(37/2,6). The calculated values of €critical for 6 « 1 are 
given in Table 7.2. 


The non-linearity of the flow law increases the range of ¢ values for which a is 
found, as Table 7.2 shows. The velocity decrease with respect to the velocity at the 


bed also gets more pronounced with increasing n (cf. Fig. 6.10). 


For « = 0, aay is situated at Z = 1. With increasing « Ut, moves downwards, 
reaches the bed and disappears for € = €critical(37/2,5). Increasing n has the effect 


min 


of moving U3" /2 Slightly further away from the bed. 
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n Ecritical (37/2, 0) 
1 1.20 
2 1.37 
3 1.45 
4 1.50 
5 1.55 


Table 7.2: €criticai(37/2,0) as a function of n for 6 < 1. All values based on numerical 
calculations. The velocity minimum above the trough of the sine curve (kx = 37/2) exists for 
€ < Ecritica The errors in the determination of Ecriticas are estimated to be less than 7%. 


Below Usr'?, the horizontal velocities increase with depth. If the explanation for the 
occurrence of extrusion flow at Storglacidren is correct, one must expect two further 
such regions above and below the riegel. This would be an interesting test of the 
appropriateness of this explanation. 


7.1.3 Flow separation 


Flow separation occurs for € > 1.8 for all values of n. It can happen for a perfectly 
lubricated bed as well as for no-slip boundary conditions. Based on the large number 
of analytical, experimental and computational demonstrations of its existence, flow 
separation for no-slip boundary conditions is known to be an universal feature of 
laminar flow in corners (Sherman, 1990, p. 265) as well as around sharp corners. It 
is now clear that flow separation must also be expected at roughness values larger 
than = 0.28 for smooth perfectly lubricated beds. 


A value of 1.8 for the slope parameter € corresponds to a roughness r of 0.28. This is 
a high roughness value but there is an example of an overdeepening in the Alps that 
has an even larger roughness value. This is the spectacular overdeepening found at 
Konkordiaplatz, Aletschglacier, Switzerland, which is about 1000 m long and 400 
m deep. Although the flow situation there is strongly three-dimensional it must be 
expected that the ice within this overdeepening hardly takes part in the overall flow 
of the glacier, if at all. The overlying ice most probably moves directly over the 
trough. Again this is a prediction that could be put to the test. 


7.1.4 Stresses 


The maximum of the effective stresses along the bed-line is found for ¢ + 0 at the 
inflection points of the sinusoidal curve. With increasing ¢ these maximum points 
move towards the peaks of the sine curve, and become stronger and more localized. 
A first order perturbation theory gives a distribution of the effective stresses that is 
independent of z. For a second order theory this is no longer true. 


91 i 


7.2 Sliding velocity 


The sliding velocity as a function of n, €, and 6 has been calculated. Eq. (3.23) is 
correct. For finite values of e the sliding velocity can be determined by the use of 
the Taylor coefficients listed in Table 6.1. 


A very interesting result is the fact that In u, depends linearly on n for all values of 
e and 6.! 


Furthermore the numerical calculations indicate that s depends linearly on 6. As n 
and « get larger the linear slope of the s(6) curve increases approximately according 
to Eq. (6.14). The sliding velocity is thus dependent on the glacier thickness and 
measurements of changes in sliding velocities with time for glaciers lying on hard 
beds could thus be caused by thickness changes. 


It is, however, difficult to see how these theoretical results could be tested or even 
used directly to estimate sliding velocities of glaciers. The model is in this respect 
definitely too idealized. On the other hand it is difficult to see how progress towards 
a theoretically-derived sliding law can be made without knowing the sliding velocity 
for this idealized model. Further work along similar lines is needed. 


1To calculate s as a function of n for a different value of n than that listed, one should use 
Table 6.1 to calculated the sliding function for several n values and then interpolate or extrapolate 
the results assuming In s = 4 + bn for some 4 and b. 
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APPENDIX A 


Second-order Solutions for Perfect 
Sliding With Regelation 


The following equations give the velocity components, the strain rates, and the 
pressure distribution of a highly viscous medium sliding over a sinusoidal bed in the 
presence of regelation. They were derived directly from Morland (1976a). Morland 
was mainly interested in the velocity and pressure distribution along the bed, and 
therefore only gave the corresponding expressions for z = 


The zx axis is parallel to the mean slope of the bed making an angle @ with respect 
to the horizontal. The coordinates of the bed as a function of x are denoted by 
zo(x). The bed is defined by 29 = asinkz. 


Regelation is only important at wavelengths comparable to or smaller than the 
transition wavelength A,, where 


2 r 7 
i= ea (A.1) 


Using the numerical values for the physical parameters from Table A.2 one finds 
that the transition wavelength is of the order of 0.5m. For a non-linear medium 
with n > 1 the transition wavelength becomes smaller (Kamb, 1970; Fowler, 1979; 
Fowler, 1981). Table A.1 compares the notation that is used here with the notation 
of several other authors. 


There are two parameters that describe the relative importance of regelation to 
viscous flow. One of them is denoted by w and is the ratio of the bed wavelength 


this thesis Nye (1969,1970) Kamb (1970) Lliboutry (1987b) Morland (1976a) 


ky k, lo W. 1/d. 
As 2n/k. Xo 2n/w, 2m, 
L L A pL L 
k k h w k 
r Vv2r ¢ e/2a 





Table A.1: Notation used in this thesis and that used by several different authors. 
k, is the controlling wavenumber, and \, the transition wavelength, A, = 2r/k,, with ky = 


aero r = a/X is the single wavelength roughness, and ¢ = ak = 2rr is the (local 
nCgn CAG \p . : 
bed) slope number. L is the latent heat of fusion per unit volume of ice. 
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Co 0.8 x 107-7 °Cm?N7! 
Q 4x 10°? Jm-?s"! 
L 2.8 x 10° J m=? 

n 3x 10!?Nm7?s 
K; 2.0Jm-!s-!°Cuh 

Kg ~~] K, 


Table A.2: Numerical values of physical parameters, based on Paterson (1981) and Hutter 
(1983). Cp is the Clausius-Clapeyron constant, Q the geothermal heat flux, L the latent heat 
of fusion per unit volume of ice, 7; is the (linear) viscosity of ice, and K, and Kg are the 
thermal conductivities of the ice and bed respectively. 


to the transition wavelength, :.e. 


w= AYA)? Sk / ky (A.2) 
L 
2k2nCo (K, + Kg) 


k, is the controlling wavenumber, given by 


9 L 


k= ——————. A.4 
2nCo (K; + Kp) oe 


Since the ratio w?/(w? + 1) will turn up in most of the following equations it will 


be abbreviated by {), or 
3? 


A aT 
This also agrees with Nye’s notation. In the no-regelation limit @; = 1, and in the 
pure-regelation limit {9 = 0. Inw is a convenient measure of the importance of 


regelation. 


(A.5) 


The second parameter is Morland’s A parameter, which will here be denoted by Ayn. 


A, = (Kit Ka) cos(a)pg Co +2@ (A.6) 
Lu 
_ akin (UR Raorle) 22), (A.7) 
L h Tp 





The effect of freezing and melting on the flow field is negligible if A,, <1, which is 
almost always the case (Morland, 1976a). 


The basal sliding velocity uy is 





T (Ww? +1) 
Tb 1 1 k 
os a eels A.9 
ne? (E+E z| (A.9) 


94 


In the no-regelation limit k/k, < 1 and 


Tp 


In the pure-regelation limit k/k, >> 1 and 
™k 
Us = eek? (A.11) 


up, has a minimum at k = k, and the largest part of the drag is contributed by 
the Fourier components of the bed with wavelengths (having the same amplitude to 
wavelength ratio) around A, (Nye, 1969). 
Note that 

i= (A.12) 

b KI nerk’ 5 

which is a useful relation that can be used to eliminate the sliding velocity from the 
following equations. 


The horizontal velocity field is given by 


=o Th 
Vz(Z,z2) = Upt on (1 — (1 — z/h)?) (A.13) 


+ up, kze* (sin(kx) — Ap, cos(kx)) € 
+ uy be~?* (cos(2kxr) + Am sin(2 kx)) (1/4 — kz/2) &? 
+ O(e?). 

Note also that for z = 0 the first order term vanishes. 


The vertical velocity field is given by 


v,(z,z) = up, f,e* (cos(kz) + Am sin(kxr)) (1+ kz) € (A.14) 
+ t 3, kze~?** (sin(2 kx) — A,, cos(2 kz)) €? 
+ O(e%). 


The corresponding strain rates are given by 


xe = Up f, zk?e~** (cos(kx) + Am sin(kz)) € (A.15) 
+ up By ke7?* (sin(2kxr) — Ap cos(2kx)) (kz — 1/2) €? 
+ O(e%). 

éxz = ; (1 — z/h) a?k* uy (A.16) 
+ uy QB, zk2e7* (Am cos(kr) — sin(kx)) € 
+ up 3, ke7?*™ (kz — 1/2) (cos(2 kx) + Am sin(2 kr)) Ee 
+ O(e%). 
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The second invariant of the strain-rate tensor is given by 


1 


te gar ne ass (A.17) 
+ u? B, (1 —2/h) (Am cos(kx) — sin(kx)) zk°a?e7** € 
so 42 Be (ie 1) ke? ** ¢ 2 
+ u2B, (1 —2z/h) (kz — 1/2) (Am sin(2 kx) + cos(2 kz)) a2k4e7?* 6? 
+ ub 6? (A? +1) (2kz — 1)sin(ke) zk8e73** 63 
+ up BP (kz — 1/2)” (4%, ea 1) pee 4ke 4 
+ Of). 


The pressure is given by 


P(z,z) = Pa + pg cos(a) (1 — 2/h) (A.18) 
us 3, nke~*? (cos(kx) + Am sin(kz)) € 

us B, nke-?*? (sin(2 kz) — Am cos(2kzx)) €? 

O(e?). 


+++ 


These expressions can be used to calculate the flow and the sliding velocity for a 
general bed geometry. Different Fourier components can simply be added linearly 
together. Calculating the regelation velocity and the viscous flow velocity separately 
for a general bed geometry and afterwards adding the two components together 
would on the other hand give an incorrect result. The reason for this is, as Nye (1969) 
explains, that the pressure distribution that results from the process of regelation 
is in general different from the one that is caused by viscous flow over a bump, so 
that the two flows will interfere. 
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